Methods and systems for genome-scale kinetic modeling

ABSTRACT

Embodiments of the present invention generally relate to the construction, analysis, and characterization of dynamical states of biological networks at the cellular level. Methods are provided for analyzing the dynamical states by constructing matrices using high-throughput data types, such as fluxomic, metabolomic, and proteomic data. Some embodiments relate to an individual, while others relate to a plurality of individuals.

RELATED APPLICATIONS

The present application is a continuation of U.S. application Ser. No.12/918,317, filed on Aug. 18, 2010, which is a national stageapplication under 35 U.S.C. §371 of PCT Application No.PCT/US2009/034587, filed on Feb. 19, 2009, which published in English asWO 2009/105591 A3 on Aug. 27, 2009, and which claims benefit of U.S.Provisional Application No. 61/029,891, entitled Methods and Systems forGenome-scale Kinetic Modeling, filed Feb. 19, 2008, the entire contentsof which applications and publication are herein incorporated byreference in their entirety.

STATEMENT OF GOVERNMENT SUPPORT

This invention was made with government support under grant numberBNGBP47 awarded by NIH. The government has certain rights in theinvention.

BACKGROUND OF THE INVENTION

1. Field of the Invention

Embodiments of the present invention generally relate to theconstruction, analysis, and characterization of dynamical states ofbiological networks at the cellular level. Methods are provided foranalyzing the dynamical states by constructing matrices usinghigh-throughput data types, such as fluxomic, metabolomic, and proteomicdata. Some embodiments relate to an individual, while others relate to aplurality of individuals.

2. Description of the Related Art

There have been significant advances in the development of statisticalanalysis of genome-scale networks (Slonim, 2002), that have beenpropelled by the availability of genome-scale high-throughput datasetsand the successes of constraint-based modeling approaches (Kummel etal., 2006; Palsson, 2006; Pharkya et al., 2004; Price et al., 2004; Reedet al., 2006). The foundation of such genome-scale analysis is builtupon the descriptions of the biochemical transformations in a network ina self-consistent and chemically accurate matrix format. Much progresshas been made with the genome-scale network reconstruction process, anda growing number of genome-scale metabolic reconstructions are nowavailable (Feist et al., 2007; Jamshidi and Palsson, 2007; Oh et al.,2007; Reed et al., 2006; Resendis-Antonio et al., 2007).

Reconstructions of genome-scale biochemical reaction networks (Reed etal., 2006) have been analyzed by topological (Barabasi and Oltvai, 2004)and constraints-based (Price et al., 2004) methods, but dynamic models,at this scale, have been unachievable to date. It turns out that thematrix representation of the biochemical transformations in a networkare not only a requisite for dynamic models but also a major determinantin their properties, and thus, it is important to have well curatedreconstructions available. The growing availability of metabolomic andfluxomic data sets (Breitling et al., 2007; Goodacre et al., 2004;Hollywood et al., 2006; Sauer, 2006) and methods to estimate thethermodynamic properties (Henry et al., 2007; Henry et al., 2006;Mavrovouniotis, 1991) of biochemical reactions has opened up thepossibility to formulate large-scale kinetic models.

The growing amount of disparate data sets (including proteomic,fluxomic, and metabolomic) has produced the unmet need to integrate andanalyze this data in a functional context that is tailored to specificindividuals. The availability of annotated genomes enabled thereconstruction of genome-scale networks, and now the availability ofhigh-throughput metabolomic and fluxomic data along with thermodynamicinformation opens the possibility to build genome-scale kinetic models.This invention describes a framework for building and analyzing suchmodels, thus satisfying those needs.

SUMMARY OF THE INVENTION

The structure of the workflow that leads to large-scale dynamic modelsis described in this invention along with the description of theframework which enables the integration of data to constructgenome-scale kinetic models.

A description of the types of data matrices which are needed toconstruct large-scale kinetic models of biological processes isprovided.

A description of the way to construct the data matrices needed for theformulation of kinetic models is provided. A description of importantproperties of these data matrices which enable approximations of thevalues in the data matrices when partial or incomplete data is availableis also provided.

A method for the analysis of kinetic models and for providing a methodfor interpreting and predicting physiological and pathophysiologicalstates of biological reaction networks is provided.

A method is provided for developing data-driven dynamic models ofbiological networks, the method comprising providing a data structurerelating a plurality of biological network reactants to a plurality ofbiological network reactions, wherein each of the biological networkreactions comprises one or more biological network reactants identifiedas a substrate of the reaction, one or more biological network reactantsidentified as a product of the biological network reaction and astoichiometric coefficient relating the substrate and the product;providing a constraint set for a plurality of biological networkreactions; providing a constraint set for the plurality ofconcentrations associated with the plurality of biological networkreactions; providing a constraint set for a plurality of kineticconstants associated with the plurality of biological network reactions;providing a constraint set for a plurality of equilibrium constantsassociated with the plurality of biological network reactions;integrating the sets of constraints, thereby forming the stoichiometric,gradient, and Jacobian matrices; providing an objective function; andpredicting a physiological function related to a biological network bydetermining at least one flux distribution that minimizes or maximizesthe objective function when the constraint set is applied to the datastructure.

A method is provided wherein the constraints for the plurality ofreactions, the constraints for the plurality of concentrations, theconstraints for the plurality of kinetic constants, and the constraintsfor the plurality of equilibrium constants are uniquely defined.

A method is provided wherein providing a constraint set for a pluralityof concentrations comprises measuring metabolite levels at steady state.

A method is provided, wherein providing a constraint set for a pluralityof kinetic constants comprises solving a steady state mass conservationrelationship comprising a stoichiometric matrix and flux vector.

A method is provided further comprising defining a unique set ofmatrices using the objective function, wherein the constraints within atleast one of the sets selected from the reactions constraint set, theconcentration constraint set, the kinetic constant constraint set, andthe equilibrium constraint set are not unique across the set.

A method is provided further comprising defining a unique set ofgradient and Jacobian matrices using order of magnitude approximations,wherein the constraints within at least one of the sets selected fromthe reactions constraint set, the concentration constraint set, thekinetic constant constraint set, and the equilibrium constraint set arenot unique across the set.

A method is provided further comprising determining a set of gradientand Jacobian matrices corresponding to the plurality of constraintconditions, wherein the constraints within at least one of the setsselected from the reactions constraint set, the concentration constraintset, the kinetic constant constraint set, and the equilibrium constraintset are not unique across the set.

A method is provided wherein constraints within at least one of theconcentration constraint set and the kinetic constant constraint set arenot unique across the set and are estimated by order of magnitudeapproximations.

A method is provided wherein the stoichiometric, gradient, and Jacobianmatrices refer to a biological network process.

A method is provided wherein the biological network process comprises atleast one of metabolism, signal transduction or signaling,transcription, and translation.

A method is provided wherein the sets of constraints correspond to anindividual.

A method is provided wherein the sets of constraints correspond to agroup or plurality of individuals.

A method is provided wherein the model is used to determine the effectsof changes from aerobic to anaerobic conditions.

A method is provided wherein the biological network reactions areobtained from a metabolic reaction database that includes thesubstrates, products, and stoichiometry of a plurality of biologicalreactions.

A method is provided for the analysis of disease states due to geneticand or environmental influences, the method comprising providing a datastructure relating a plurality of reactants to a plurality of reactions,wherein each of the reactions comprises one or more reactants identifiedas a substrate of the reaction, one or more reactants identified as aproduct of the reaction and a stoichiometric coefficient relating thesubstrate and the product; providing a constraint set for the pluralityof reactions; providing a constraint set for a plurality ofconcentrations associated with the plurality of reactions; providing aconstraint set for a plurality of kinetic constants associated with theplurality of reactions; providing a constraint set for a plurality ofequilibrium constants associated with the plurality of reactions;integrating the sets of constraints, thereby forming the stoichiometric,gradient, and Jacobian matrices; providing an objective function;modifying at least one of the reaction constraint set, the concentrationconstraint set, the constraint set for a plurality of kinetic constantsand the set for a plurality of equilibrium constants based at leastpartly on altered environmental conditions; and predicting aphysiological function related to the disease state by determining atleast one flux distribution that minimizes or maximizes the objectivefunction when at least one constraint set is applied to the datastructure.

A method is provided wherein the model is used to determine the dynamicsand consequences of genetic defects selected from the group consistingof deficiencies in metabolic enzymes and metabolic transporters.

A method is provided wherein the genetic defects comprise thefunctioning of a compound selected from the group consisting ofphosphofructokinase, phosphoglycerate kinase, phosphoglycerate mutase,lactate dehydrogenase adenosine deaminase, ABC transporters, the SLCclass of transporters, and the cytochrome P450 class of enzymes.

A method is provided wherein providing a constraint set for a pluralityof concentrations comprises measuring metabolite levels at steady state.

A method is provided wherein providing a constraint set for a pluralityof kinetic constants comprises solving a steady state mass conservationrelationship comprising a stoichiometric matrix and flux vector.

A method is provided wherein the modified sets of constraints correspondto an individual.

A method is provided wherein the modified sets of constraints correspondto a group or plurality of individuals.

A method is provided wherein a therapeutic regime is applied to regulatea physiological function based on at least one of the modified setsbased at least partly on genetic variations.

A method is provided wherein a therapeutic regime is applied to regulatea physiological function based on at least one of the modified setsbased at least partly on altered environmental conditions.

A method is provided for the analysis of therapies for drug treatments,the method comprising providing a data structure relating a plurality ofreactants to a plurality of reactions, wherein each of the reactionscomprises one or more reactants identified as a substrate of thereaction, one or more reactants identified as a product of the reactionand a stoichiometric coefficient relating the substrate and the product;providing a constraint set for the plurality of reactions; providing aconstraint set for a plurality of concentrations associated with theplurality of reactions; providing a constraint set for a plurality ofkinetic constants associated with the plurality of reactions; providinga constraint set for a plurality of equilibrium constants associatedwith the plurality of reactions; integrating the sets of constraints,thereby forming the stoichiometric, gradient, and Jacobian matrices;providing an objective function; providing a experimental constraint setfor the plurality of reactions, concentrations, or kinetic constants foran altered set of constraints due to the treatment with a drug orchemical agent; and predicting a physiological function related to adrug treatment by determining at least one flux distribution thatminimizes or maximizes the objective function when a constraint set isapplied to the data structure.

A method is provided wherein a set of fluxes and concentrations aremeasured in an individual before and after giving a dose of atherapeutic agent to the individual and constructing kinetic networkmodels according to the measured data specific to the individual.

A method is provided wherein providing a constraint set for a pluralityof concentrations comprises measuring metabolite levels at steady state.

A method is provided wherein providing a constraint set for a pluralityof kinetic constants comprises solving a steady state mass conservationrelationship comprising a stoichiometric matrix and flux vector.

A method is provided wherein a physiological function related to a drugtreatment is instead predicted by determining at least one concentrationthat maximizes or minimizes the objective function when the experimentalconstraint set is applied to the data structure.

A method is provided wherein a drug or chemical agent is applied toregulate a physiological function based on the at least one experimentalconstraint set for an altered state of constraints due to the treatmentwith a drug or chemical agent.

A method is provided wherein a drug or chemical agent is applied toregulate a physiological function based on the at least one experimentalconstraint set which maximizes or minimizes the objective function whenthe experimental constraint set is applied to the data structure.

A method is provided wherein the sets of constraints correspond to anindividual.

A method is provided wherein the sets of constraints correspond to agroup or plurality of individuals.

A method is provided for accounting for small metabolite regulation, themethod comprising: providing a data structure relating a plurality ofreactants to a plurality of reactions, wherein each of the reactionscomprises one or more reactants identified as a substrate of thereaction, one or more reactants identified as a product of the reactionand a stoichiometric coefficient relating the substrate and the product,in which the reactions may reflect regulatory processes; providing aconstraint set for the plurality of reactions; providing a constraintset for a plurality of concentrations associated with the plurality ofreactions; providing a constraint set for a plurality of kineticconstants associated with the plurality of reactions; providing aconstraint set for a plurality of equilibrium constants associated withthe plurality of reactions; integrating the sets of constraints, therebyforming the stoichiometric, gradient, and Jacobian matrices; providingan objective function; and predicting a physiological function relatedto the small metabolite regulation by determining at least one fluxdistribution that minimizes or maximizes the objective function when theconstraint set is applied to the data structure.

A method is provided wherein providing a constraint set for a pluralityof concentrations comprises measuring metabolite levels at steady state.

A method is provided wherein a physiological function related to thesmall metabolite regulation is instead predicted by determining at leastone concentration that minimizes or maximizes the objective functionwhen the constraint set is applied to the data structure.

A method is provided wherein reactants and products include proteinsthat interact with a plurality of other components.

A method is provided wherein reactants or products activate a particularenzyme or plurality of enzymes.

A method is provided wherein reactants or products inhibit a particularenzyme or plurality of enzymes.

A method is provided for determining one or more parameters tocharacterize different functional states, the method comprising:analyzing a network function; determining time scales of interest; andidentifying parameters that describe the network function.

A method is provided wherein the parameters comprise at least one of areaction, a concentration and a rate constant.

A method is provided further comprising characterizing the dynamicbehavior of the network using the time scale determination; andidentifying reactions or concentrations that determine the networkfunction.

A computer readable medium is provided, the medium containing softwarethat, when executed, causes the computer to perform the acts of:providing a data structure relating a plurality of reactants to aplurality of reactions, wherein each of the reactions comprises one ormore reactants identified as a substrate of the reaction, one or morereactants identified as a product of the reaction and a stoichiometriccoefficient relating the substrate and the product; providing aconstraint set for the plurality of reactions; providing a constraintset for a plurality of concentrations associated with the plurality ofreactions; providing a constraint set for a plurality of kineticconstants associated with the plurality of reactions; providing aconstraint set for a plurality of equilibrium constants associated withthe plurality of reactions; integrating the sets of constraints, therebyforming the stoichiometric, gradient, and Jacobian matrices; providingan objective function; and predicting a physiological function bydetermining at least one flux distribution that minimizes or maximizesthe objective function when the constraint set is applied to the datastructure.

A computer implemented system is provided comprising: a computingenvironment; a storage in data communication with the computingenvironment and configured to store data structure relating a pluralityof reactants to a plurality of reactions; and a software programoperating on the computing environment and configured to: provide thedata structure relating the plurality of reactants to the plurality ofreactions, wherein each of the reactions comprises one or more reactantsidentified as a substrate of the reaction, one or more reactantsidentified as a product of the reaction and a stoichiometric coefficientrelating the substrate and the product; provide a constraint set for theplurality of reactions; provide a constraint set for a plurality ofconcentrations associated with the plurality of reactions; provide aconstraint set for a plurality of kinetic constants associated with theplurality of reactions; provide a constraint set for a plurality ofequilibrium constants associated with the plurality of reactions;integrate the sets of constraints, thereby forming the stoichiometric,gradient, and Jacobian matrices; provide an objective function; andpredict a physiological function by determining at least one fluxdistribution that minimizes or maximizes the objective function when theconstraint set is applied to the data structure.

A system is provided comprising means for providing a data structurerelating a plurality of reactants to a plurality of reactions, whereineach of the reactions comprises one or more reactants identified as asubstrate of the reaction, one or more reactants identified as a productof the reaction and a stoichiometric coefficient relating the substrateand the product; means for providing a constraint set for the pluralityof reactions; means for providing a constraint set for a plurality ofconcentrations associated with the plurality of reactions; means forproviding a constraint set for a plurality of kinetic constantsassociated with the plurality of reactions; means for providing aconstraint set for a plurality of equilibrium constants associated withthe plurality of reactions; means for integrating the sets ofconstraints, thereby forming the stoichiometric, gradient, and Jacobianmatrices; means for providing an objective function; and means forpredicting a physiological function by determining at least one fluxdistribution that minimizes or maximizes the objective function when theconstraint set is applied to the data structure.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a graphical diagram showing relationships between thefundamental data matrices describing dynamic states of biologicalnetworks: the stoichiometric matrix S and the gradient matrix G.

FIG. 2 is a graphical diagram illustrating a workflow for constructingapproximations to genome-scale kinetic networks and how the fourproperties discussed in this manuscript and highlighted in FIG. 1 enablesuch reconstructions.

FIG. 3A, represented in FIGS. 3A-1 and 3A-2, shows an example of theglycolytic pathway in Homo sapiens. A reaction map of the glycolyticpathway is shown. The decomposition of the Jacobian (J_(x)) into thestoichiometric, κ, and Γ matrices follows below (1-norm used for thefactorization of Γ). The negative transpose of the gradient matrix isshown directly below the stoichiometric matrix, demonstrating thestructural similarity.

FIG. 3B, represented in FIGS. 3B-1 and 3B-2, explicitly illustrates theprinciple of duality and time scale decomposition via the resulting datamatrices. The Jacobian duals are shown; they are related by the gradientmatrix. Hierarchical analysis can be carried out of the network in termsof metabolites or fluxes. The resultant modal matrices can be related toone another via the stoichiometric matrix. As illustrated in FIG. 4,sometimes it is convenient to think about the hierarchical structure interms of metabolites and sometimes it is more intuitive to think interms of the fluxes. The network was constructed using equilibriumconstants and concentrations from a kinetic red cell model, which hasbeen validated in the literature. Network dynamics were then describedusing mass action kinetics for a particular steady state. Abbreviations:TS: time scale in hours, HK: hexokinase, PGI: phosphoglucoisomerase,PFK: phosphofructokinase, ALD: fructose-bisphosphate aldolase, TPI:triose-phosphate isomerase, GAPD: glyceraldehyde phosphatedehydrogenase, PGK: phosphoglucokinase, PGM: phosphoglucomutase, ENO:enolase, PYK: pyruvate kinase, LDH: lactate dehydrogenase, LEX: lactatetransporter, ATPase: ATP hydrolysis (demand/utilization).

FIG. 4 shows the hierarchical reduction of glycolysis.

FIG. 5 shows the hierarchical analysis and simplification of red cellmetabolism.

FIG. 6 is a block diagram of the general diagnostic and therapeuticmethodologies implementing one embodiment of the invention.

FIG. 7 is a generalized flow diagram of one embodiment of the invention.

FIG. 8 is a block diagram of one embodiment for preparing the modeldescribed herein.

FIG. 9 shows the hierarchical reduction of the pathway discussed inExample 5 below. The reaction scheme for hexokinase as described by (23)hkE represents the unbound, free form of the enzyme. The individualsteps and interactions accounting for the catalytic steps in PFK inaddition to the allosteric interactions. Reaction scheme for theRapoport-Leubering shunt which is carried out by the same enzyme (31).The reaction schemes for Glucose-6-Phosphate Dehydrogenase as describedby (23) g6pdER represents the active enzyme. The reaction schemes forAdenosine Kinase as described by (25) akE represents the enzyme in thecatalytic state and akET represents the tense or inactive form of theenzyme. Other abbreviations: G6P: Glucose-6-phosphate, F6P:Fructose-6-phosphate, FDP: Fructose-1,6-bisphosphate, DHAP:Dihydroxyacetone phosphate, GAP: Glyceraldehyde-3-phosphate, DPG 13:1,3-bisphosphoglycerate, DPG23: 2,3-bisphosphoglycerate, PG3:3-phosphoglycerate, PG2: 2-phosphoglycerate, PEP: Phosphoenolpyruvate,PYR: Pyruvate, LAC: Lactate, NADH: Nicotinamide adenine dinucleotide(reduced), GL6P: 6-Phospho-D-glucono-1,5-lactone, GO6P:6-Phospho-D-gluconate, NADPH: Nicotinamide adenine dinucleotidephosphate (reduced), GSH: glutathione (reduced), RU5P:Ribulose-5-phosphate, RSP: Ribose-5-phosphate, X5P:Xylulose-5-phosphate, S7P: Sedoheptulose-7-phosphate, E4P:Erythrose-4-phosphate, ADO: Adenosine, AMP: Adenosine monophosphate,ADP: Adenosine diphosphate, ATP: Adenosine triphohsphate, PRPP:5-Phospho-D-ribose-I-diphosphate, IMP: Inosine monophosphate, INa:Inosine, HX: Hypoxanthine, RIP: Ribose-1-phosphate, ADE: Adenine, Phosi:Monophosphate.

FIG. 10 is a tiled phase plane diagram for selected grouped variablesand reaction fluxes for the cell in response to a transient,instantaneous decrease in ATP by 60%. HK Ratio, PFK Ratio, DPGM Ratio,G6PD Ratio, and AK Ratio are the ratios of the free to bound forms ofthe corresponding enzymes. The lower left triangle depicts thetrajectories of the corresponding variables. The upper right triangleshows the correlation coefficient (when the correlation coefficient isgreater than 0.85, the slope of the line appears below it). The G6PDratio, NADPHINADP ratio, and G6PD I flux are all highly correlated, asare the ATPase flux and the ATP/ADP ratio.

FIG. 11, represented in FIGS. 11A and 11B, is a tiled phase planediagram for all of the fluxes involving the enzyme PFK for the cell inresponse to a transient, instantaneous decrease in ATP by 60%. Note thatnot all of the steps are correlated with one another. Indeed, the netflux through the enzyme is actually a sub-network in and of itself. Thecatalytic step is PFK3. PFK14 through PFK17 reflect ATP binding to thetense form of the enzyme and PFK18 through PFK21 represent Mg binding tothe tense form of the enzyme. Correlations among some of the fluxes areobserved, for example, the sequential binding of multiple ATP moleculesto the enzyme. Other fluxes move dynamically independently from theother PFK fluxes, for example PFK15 and PFKI8.

FIG. 12 is a tiled phase plane diagram for selected grouped variablesand reaction fluxes for the cell in response to a transient,instantaneous decrease in NADPH by 80%. Similar to the simulation withATP reduction, the G6PD1 flux, G6PD Ratio, and NADPHINADP ratio are allhighly correlated during the time course.

FIG. 13 is a table of the stoichiometric matrix S for the mass actionstudy of the red blood cell, which for ease of viewing is represented inFIGS. 13-1, 13-2, 13-3 and 13-4.

FIG. 14 is a table of the gradient matrix G for the mass action study ofthe red blood cell, which for ease of viewing is represented in FIGS.14-1, 14-2, 14-3, 14-4, 14-5 and 14-6.

FIG. 15 is a table of the stoichiometric matrix S for the regulatedstudy of the red blood cell, which for ease of viewing is represented inFIGS. 15-1 through 15-27.

FIG. 16 is a table of the gradient matrix G for the regulated study ofthe red blood cell, which for ease of viewing is represented in FIGS.16-1 through 16-27.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

Methods for constructing and analyzing dynamic in silico models ofbiological networks are provided herein. A dynamic model can be used tosimulate different aspects of cellular behavior of human cells underdifferent environmental and genetic conditions, thereby providingvaluable information for a range of industrial and researchapplications. Models and methods disclosed herein can also be extendedto simulate the activity of multiple interacting cells, includingorgans, physiological systems, and whole body metabolism. Developingkinetic models of biological networks can be used to inform and guidethe research process, potentially leading to the identification of newenzymes, medicines or metabolites of commercial importance.

As used herein, the term “concentration” refers to a numerical valuewith physical units of mass*volume⁻¹, such as molar and millimolar. Suchquantities include but are not limited to molecules in a biochemicalnetwork, such as glucose, pyruvate, lactate, enzymes that carry outbiochemical transformations or transport reactions such as hexokinase,adenosine deaminase, and sodium-potassium ATPase pump, or the number ofcells in a volume of tissue.

As used herein, the term “initial conditions” refers to the set ofinitial variable values that need to be specified for carrying outdynamic simulations in a system of dynamic differential equations.

As used herein, the term “dynamic mass balance” refers to a linear firstorder differential equation with respect to time which accounts for thesummed reactions producing or consuming each component in a specifiedbiological network.

As used herein the term “flux” or “reaction flux” refers to a singleflux in a flux distribution. Individual fluxes can be represented as thedifference between a forward and reverse flux, with units mass, ornumber of moles, or number of molecules per unit volume per unit time.

As used herein, the term “flux distribution” refers to a directional,quantitative list of values corresponding to the set of reactions in anetwork, representing the mass flow per unit time for each reaction inthe network being analyzed.

As used herein, the term “concentration distribution” refers to anon-negative, quantitative list of values corresponding to the set ofconcentrations of variables, including but not limited to metabolitesand enzymes, in a network, representing the mass, or number of moles, ornumber of molecules per unit volume for each variable.

As used herein, the term “biological network” refers to assembledreactions reflecting biological processes. Examples of biologicalnetworks include but are not limited to metabolic networks, signalingnetworks, and regulatory networks.

As used herein, the term “stoichiometric matrix” or “S” refers to amatrix with the stoichiometric coefficients for reactions represented bythe columns and the substrates and products in the rows. Thestoichiometric matrix may be written for any biological network.Examples include but are not limited to metabolic networks and signalingnetworks.

As used herein, the term “gradient matrix” or “G” refers to a matrixwhose row entries are the first order partial derivatives of reactionfluxes with respect to the metabolites (see equation 2). The gradientmatrix is specific for a particular biological network. The number ofrows in the gradient matrix is equal to the number of reactions in thenetwork. The number of columns in the gradient matrix is equal to thenumber of metabolites or compounds in the network.

As used herein, the term “Jacobian”, “Jacobi matrix”, “Jacobian matrix”,“J”, or “J_(x)” refers to a matrix of the first order partialderivatives of the linearized dynamic mass balances, in which thevariables are the metabolites or compounds. The Jacobian matrix isderived when a system of differential equations is linearized around areference state such that only the first order terms remain. TheJacobian is a square matrix and has the same number of rows and columnsas there are metabolites or compounds in the network. “J” can be writtenas the matrix product of S and G (see equations 1 and 3).

As used herein, the term “flux Jacobian” or “J_(v)” refers to a matrixof the first order partial derivatives of the linearized dynamic massbalances, in which the variables are the fluxes. The flux Jacobian is asquare matrix and has the same number of rows and columns as there arereactions or fluxes in the network. The flux Jacobian can be written asthe matrix product of G and S (see equation 5).

As used herein, the term “kappa”, “kappa matrix”, or κ refers to adiagonal matrix composed of kinetic constants. This matrix can becalculated in a variety of ways. For example the rows of G can be scaledor normalized. The factorization may be carried out but not limited tothe following manner, by factoring out the 1-norm, 2-norm, or sum of therate constants with non-zero entries in the row.

As used herein, the term “gamma”, “gamma matrix”, or Γ refers toresulting matrix after κ has been factored out of G.

As used herein, the term “transpose”, transpose of a matrix, or^(T)refers to a standard mathematical operation in which the elements of amatrix are reflected about the diagonal entries of the matrix, such thatfollowing the transformation the entries in the rows of the matrix nowappear in the columns of the matrix and the entries in the columns ofthe matrix appear in the rows of the matrix.

As used herein, the term “time constant” refers to a numerical valuewith physical units of time, such as days, hours, minutes, andmicroseconds. Examples of time constants include but are not limited tothe negative reciprocal of eigenvalues or the reciprocal entries of thekappa matrix, which may correspond to a particular reaction or set ofreactions.

As used herein, the term “kinetic constant” or reaction rate constantrefers to a time constant characteristic of a particular reaction or setof reactions. The physical units are [concentration]^(n+1)*time⁻¹, wheren is the kinetic order of the corresponding reaction. For example, azero order, first order, and second order kinetic constant has units[concentration]*time⁻¹, time⁻¹, and concentration⁻¹*time⁻¹ respectively.

As used herein, the term “steady state” refers to any of a number ofconditions for which the flux distribution and concentrationdistribution do not change over time. Equilibrium is a special casewhich also satisfies the steady state. In general, unless otherwisespecified, references to the steady state imply a non-equilibrium steadystate.

As used herein, the term “concentration deviation variable”, “deviationvariable”, or x′ is vector equal to the concentration variable vectorminus a reference state. This reference state may be one of any numberof concentration distributions including but not limited to steady stateconcentrations.

As used herein, the term “time scale decomposition” refers to theestablished practice of characterizing network dynamics by decomposingthe Jacobian matrix into its constituent eigenvectors and eigenvalues,from which dynamically independent modes and time constants can bederived.

As used herein, the term “modal matrix” or “modes” refers to the inverseof the eigenvector matrix of J.

The reconstruction of a genome-scale reaction network requires theidentification of all its chemical components and the chemicaltransformations that they participate in. This process primarily relieson annotated genomes and detailed bibliomic assessment. See Reed et al.,2006, which is hereby incorporated by reference in its entirety. Theresult of this process is the stoichiometric matrix, S, that is used inthe dynamic mass balances

dx/dt=S·v(x)x=x ₀ at time t=0  (1)

that are the basis of all kinetic models. Here d(.)/dt denotes the timederivative, x is the vector of the concentrations of the compounds inthe network and v(x) is the vector of the reaction rates. Allbiochemical transformations are fundamentally uni- or bi-molecular. Suchreactions can be represented by mass action kinetics, or generalizationsthereof. See Segel, 1975, which is hereby incorporated by reference inits entirety. The net reaction rate for every elementary reaction in anetwork can be represented by the difference between a forward andreverse flux (e.g. see FIG. 1). This commonly used formulation is basedon several well-known assumptions, such as constant temperature, volume,and homogeneity of the medium. If S, v(x) and the initial conditions(x₀) are known, then these ordinary differential equations can benumerically solved for a set of conditions of interest.

FIG. 1 shows a structural similarity between S and G^(T), highlighted bythe fact that their row reduced forms are equivalent and illustrated fora simple example involving two reactions. The correspondingstoichiometric matrix and gradient matrix (for mass action kinetics) isshown. The decomposition of the Jacobian matrix into the stoichiometricand gradient matrices is biologically meaningful and driven by theunderlying chemistry, kinetics, and thermodynamics. There are differentways to factor G into Γ and the diagonal κ matrix. In the embodimentshown in FIG. 1, the forward rate constant is factored out of G. Aduality exists between the concentrations and fluxes. The practicalsignificance of this complementary relationship is that in the linearregime, the relationship between fluxes and concentrations can bedetermined independently of the specific rate law formulation, if onecan approximate the gradient matrix. When certain reactions occur muchfaster than others, the related metabolites pool together into aggregatevariables. For the example in A, when the forward and reverse rateconstants of the first reaction (2A⇄B; dimerization reaction) are muchfaster than the rate constants for the second reactions (B+X⇄C+Y;cofactor exchange reaction), then the substrate and products of thefirst reaction form an aggregate pool. Since B is a ‘dimer’ of A, thepool will consist of A and B in a 1:2 ratio. Furthermore, since thefirst reaction occurs much more quickly than the second, the interactionbetween the A+2B pool and the C(X and Y are considered cofactors in thisexample) is determined by the rate constant(s) for the second reaction.

The characterization of the dynamic states of networks can be studiedthrough numerical simulation or through using mathematical analysis. Asimulation is context-dependent and represents a case study.Mathematical methods for the analysis of model characteristics typicallyrely on studying the properties of the transformation between theconcentrations and fluxes. The analysis of such fundamental propertiesnormally relies on the linearization of the governing equations at adefined condition. The linearization of the dynamic mass balanceequations comes down to the linearization of the reaction rate vector toform the gradient matrix

G=[g _(ij) ]=[∂v _(i) /∂x _(j)]  (2)

in which ∂·/∂x_(j) is the partial derivative operator. The Jacobianmatrix at a reference state X_(ref) can then be described:

dx′/dt=S·G·x′=J·x′  (3)

where x′=x−x_(ref) and J is the well-known Jacobian matrix. Analysis ofthe characteristics of Jacobian matrix is standard procedure inmathematical analysis of system dynamics (e.g. Strogatz, 1994, which ishereby incorporated by reference in its entirety). The application ofthese methods to biochemical networks has been carried out for decades.See Heinrich et al., 1977; Heinrich et al., 1991; Heinrich and Sonntag,1982; Palsson et al., 1987, all of which are hereby incorporated byreference in their entireties. In recent years there has been renewedinterest and recently further developments in the dynamic analysis ofthe properties of J have appeared. See Bruggeman et al., 2006; Famili etal., 2005; Grimbs et al., 2007; Kauffman et al., 2002; Steuer et al.,2006; Teusink et al., 2000, all of which are hereby incorporated byreference in their entireties.

The integration of genomic, bibliomic, fluxomic, and metabolomic datawith thermodynamic information into dynamic models of metabolism isillustrated in FIG. 2. The process of reconstructing S is now welldeveloped. See Palsson, 2006, which is hereby incorporated by referencein its entirety. The formulation of G can now be performed based onmetabolomic data and methods to estimate thermodynamic properties. Thechemical equations that make up S determine the location of the non-zeroelements in G and hence its structure.

The availability of genomic and bibliomic data can be used to define S(FIG. 1), which has been described in detail (Palsson, 2006). Theincreasing availability of metabolomic data and approximations ofthermodynamic parameters, such as the equilibrium constants forreactions, will enable the definition of G for various networks. Thedefinition of both S and G subsequently enable the calculation of J andJ_(v). The definition of the concentration and flux Jacobian matricesprovide the basis for the characterization and analysis of the dynamicproperties of a network.

Table 1A compares properties of S and G, highlighting importantfeatures. For example, the formation of pools and reaction co-sets isdetermined by S. See Heinrich and Schuster, 1996; Jamshidi and Palsson,2006; Papin et al., 2004, all of which are incorporated by reference intheir entireties. However, time scale separation is in the realm of G.This decomposition factors the various underlying data needed for theformulation of genome-scale kinetic models. Furthermore, it illustratesthe various underlying disciplinary interests that need to be consideredand integrated to successfully achieve the analysis of dynamic states ofgenome-scale networks.

TABLE 1A A comparison of the properties of the stoichiometric matrix andgradient matrix. Properties Stoichiometric Matrix Gradient Matrixbiological species differences individual differences distal causationproximal causation genetic genomic characteristics geneticcharacteristics represents a species represents an individual informaticannotated genome kinetic data bibliomic metabolomics comparativegenomics fluxomics physico-chemical chemistry kinetics conservation lawsthermodynamics mathematical integer entries real numbers knowable matrixentries have errors numerical sparse sparse well-conditionedill-conditioned non-stiff leads to stiffness systemic pool formationtime scale separation network structure dynamic function

As mentioned supra, simulation is context-dependent and the respectivevalues will depend on the domain being modeled.

For example, in the bioinformatic context, S is primarily derived froman annotated genomic sequence fortified with any direct bibliomic dataon the organisms' gene products. The construction of G will rely onfluxomic and metabolmic data, in addition to direct kineticcharacterization of individual reactions and assessment of thermodynamicproperties.

In the physico-chemical context, S represents chemistry (i.e.stoichiometry of reactions), while G represents kinetics andthermodynamics. The chemical information is relatively easy to obtain,the thermodynamics information harder but possible (see Henry et al.,2007, which is incorporated by reference in its entirety), and thekinetic information is the hardest to acquire. The former two representhard physico-chemical constraints, while the third representsbiologically manipulable numbers; i.e., reaction rates are acceleratedby enzymes.

In the biological and genetic considerations context, the matrix S isreconstructed based on the content of a genome and is a property of aspecies. It has thus been used productively for the analysis of distalcausation (Fong and Palsson, 2004; Harrison et al., 2007; Ibarra et al.,2002; Pal et al., 2006, which are all incorporated by reference in theirentireties). Distal (or ultimate) causation results from (genomic)changes that occur from generation to generation, in contrast toproximal (or proximate) causation which occurs against a fixed geneticbackground (i.e. an individual) (Mayr, 1961, which is incorporated byreference in its entirety). G is genetically derived and can representthe variations amongst individuals in a biopopulation. It is importantin studying proximal causation and the differences in phenotypes ofindividuals in a biopopulation. For example, many disease states inhigher order organisms result from disruptions or deficiencies in enzymekinetics (Jamshidi et al., 2002, which is incorporated by reference inits entirety). These changes are reflected in G since this contains theinformation about kinetics, consequently the analysis of disease states,inter-individual differences, and transitions from a healthy to diseasestate in a particular individual will generally focus on G.

Generally speaking, the equations (S) will represent the pertinentphysical laws and the matrices (G) will represent data. Whereas S is a‘perfect’ matrix comprised of integers (i.e., digital), G is an analogmatrix whose entries are real numbers, and values of G may only be knownwithin an order of magnitude. From a numerical and mathematicalstandpoint, S is a well conditioned matrix comprised of integers(−2,−1,0,1,2), whereas G is an ill-conditioned matrix of real numbersthat can differ up to 10 orders of magnitude in their numerical values.This property leads to time scale separation. Both matrices are sparse;that is most of their elements are zero.

The elements of G, in the absence of comprehensive fluxomic,metabolomic, or thermodynamic information, may be approximated to theappropriate order of magnitude in order to enable the construction of akinetic model for analysis.

The wide differences in the numerical values of the entries of G leadsto its factorization (G=κ·Γ) by scaling it by the length of its rows toyield a factorization of J into three matrices:

J=S·κ·Γ  (4)

where κ is a diagonal matrix of the norms of the rows in G (FIG. 1.B.1).We emphasize that the i_(th) column of S contains the stoichiometriccoefficients of the i_(th) reaction in the network and the i_(th) row inG contains the forward and reverse rate constants of that same reaction.Thus, the reciprocal of the diagonal entries, 1/κ_(ii), correspond tocharacteristic time constants of the corresponding reactions. Theirnumerical values will differ significantly.

The factorization of the Jacobian in equation (4) shows that the studyof the dynamic properties of biochemical networks can be formallydecomposed into chemistry (represented by S), kinetics (represented byκ), and thermodynamic driving forces (represented by Γ). The effects ofeach can thus be formally determined. Chemistry and thermodynamics arephysico-chemical properties, whereas the kinetic constants arebiologically set through a natural selection process. The particularnumerical values (chosen through the selection process) lead to theformation of biologically meaningful dynamic properties of the network.These biological features of the network can be assessed through timescale decomposition.

The glycolytic pathway (FIG. 3A) can be used to illustrate this and theproperties below. A full kinetic model of red cell metabolism (Jamshidiet al., 2001) is available and the stoichiometric and gradient matricesare readily obtained for its glycolytic pathway. G can be factored intoκ and Γ (FIG. 3A). The elements in K are spread over approximately ten(log [κ_(max)/κ_(min)]=9.7) orders of magnitude. The matrix S isuniversal and so are the elements of G for a given set ofphysico-chemical conditions, such as temperature

It follows from Equation (2) that if s_(ij)=0 then g_(ji) is also zero;that is if a compound does not participate in a reaction, it has nokinetic effect on it. Conversely, if s_(ij)≠0 then g_(ji) is also notzero. Further inspection reveals the property that S is structurallysimilar to −G^(T) as illustrated in FIG. 1.B.2. Thus, the non-zeroentries in S have corresponding non-zero elements in −G^(T), but with adifferent numerical value. This fundamental feature shows that thetopology of the network as reflected in S has a dominant effect on itsdynamic features; providing another example of the biological principlethat structure has a dominant effect on function.

The duality relationship illustrates that fluxes or concentrations canbe used as the independent variables: A flux deviation variable, v′ canbe defined such that v′=G·x′, from which it follows that

dv′/dt=G·S·v′=J _(v) ·v′  (5)

This transformation illustrates the switch from concentrations to fluxesas the independent variables. Fluxes can tie together the multiple partsof a network to form its overall functions. Furthermore, the ability torelate the fluxes and concentrations independently of a specific ratelaw formulation, if the elements of G can be approximated, hassignificant implications for the construction and analysis of kineticnetworks.

The two Jacobian matrices, G·S and S·G, are similar and shareeigenvalues. Equation (5) is the ‘dynamic’ flux balance equation sincethe variables in it are the fluxes, v′. One can thus analyze networkproperties either in terms of concentrations or fluxes as illustratedwithin FIG. 1.B.3. The fluxes are ‘network’ variables, as they tie allthe components together, while the concentrations are ‘component’variables.

The properties of the Jacobian matrix that determine the characteristicsof the network dynamics are its eigen-properties. The eigenvalues arenetwork-based time constants (in contrast to the reaction based timeconstants in κ). Formally, the standard eigen-analysis is performed bythe diagonalization of the Jacobian matrix as:

J=M·Λ·M ⁻¹  (6)

where M is the matrix of eigenvectors, Λ is a diagonal matrix of theeigenvalues, and M⁻¹ is the matrix of eigenrows. If the decomposition ofequation (6) is introduced into equation (3) we obtain differentialequations for the modes (m=M⁻¹·x′)

dm/dt=Λ·m, or dm _(i) /dt=λ _(l) ·m _(i), or m _(i)(t)=m _(i,0)exp(λ_(l)t)  (7)

or a set of completely decoupled dynamic variables; that is, each m_(i)moves on it own time scale defined by λ_(i) independent of all the otherm_(j). The eigenrows give lumped or aggregate variables that moveindependently on each time scale since m is a set of dynamicallyindependent variables. The eigenvectors form a natural coordinate systemto describe the dynamics motion of the modes. We note that thisdecomposition can be applied to J or J_(v) (FIG. 3B.4). The eigenvalueswill be the same whereas the eigenvectors and eigenrows will differsince the variable sets (concentrations vs. fluxes) will not be thesame.

The distribution of the numerical values of the eigenvalues is the basisfor time scale separation. Time scale separation forms the basis fordecomposition of biochemical reaction networks in time and theinterpretation of the biochemical events that take place on the varioustime constants. Time scale separation has been analyzed in the contextof biological networks and shown to lend insight and enable thesimplification of theses networks (Heinrich and Sonntag, 1982; Palssonet al., 1987; Reich and Selkov, 1981, which are all incorporated byreference in their entireties).

Methods described herein enable the explicit construction of kineticmodels using the incorporation of different data types. The properties,completeness, and accuracy of the data can be explicitly traced todynamic properties. This decomposition will not only tie the modelsdirectly back to the data, but also explicitly gives us the parts of themodel that are under biological control and subject to change withadaptation or evolution. Measurement uncertainties are primarily in κand are subject to evolutionary changes. These ‘biological’ designparameters will likely need to be dealt with through the use of methodsthat bracket the range of values within uncertainty limitations.

Embodiments show how data can be used for the full delineation of thechemical equations that make up a network and ideally theirrepresentation as net combinations of elementary reactions (i.e.v_(net)=v⁺−v⁻). In this format the structure of the gradient matrix canbe determined, but integration of multiple networks is also possible andthe explicit analysis of the effects of regulatory molecules is alsoenabled. Furthermore, the underlying bilinear kinetic nature of networkdynamics is explicitly recognized by embodiments disclosed herein, asall chemical reactions are combinations of bilinear interactions,including regulatory processes (Segel, 1975).

The duality between flux distributions and concentration distributionsenabled by the transformations of the data matrices, G·S and S·G allowsrepresentation of network kinetics in terms of concentration variablesor flux variables. The systems biology paradigm of‘components→networks→computational models→physiological states’naturally leads to the use of fluxes as variables to characterize thefunctional states of a network. Fluxomic data ties components in anetwork together and is interpreted through network models, whereasconcentration data is component data. Fluxes have been widely used forsteady state analysis and can now be used to study dynamic states aswell.

Time scale decomposition has been known and studied for several decades(Heinrich et al., 1977; Heinrich and Sonntag, 1982; Palsson et al.,1987; Reich and Selkov, 1981, which are all hereby incorporated byreference in their entireties). Such studies have primarily beenperformed for small-scale models today, but their conceptual foundationhas been established. At larger scales, new issues will arise. These arelikely to include, data sensitivity, course-graining and modularizationof physiological functions in time. New methods to study the basesvectors in M and M⁻¹ that directly relate them to biochemistry andphysiological functions need to be established. The promise of theelucidation of (dynamic) structure-(physiological) functionrelationships (Palsson et al., 1987, which is hereby incorporated byreference in its entirety) may now grow to large-scale networks.

As an example, models derived by methods disclosed herein can be used todetermine the effects of changes from aerobic to anaerobic conditions,such as occurs in skeletal muscles during exercise or in tumors, or todetermine the effect of various dietary changes. These modes can also beused to determine the dynamics and consequences of genetic defects, suchas deficiencies in metabolic enzymes such as phosphofructokinase,phosphoglycerate kinase, phosphoglycerate mutase, lactate dehydrogenaseand adenosine deaminase. As another example, the dynamics andconsequences of defects in metabolic transporters, such as ABCtransporters, the SLC class of transporters, as well as the cytochromeP450 class of enzymes, may be determined.

The kinetic network models can also be used to choose appropriatetargets for drug design in individuals or groups of individuals. Suchtargets include genes, proteins or reactants, which when modulatedpositively or negatively in a simulation produce a desired therapeuticresult. Models and methods disclosed herein can also be used to predictthe effects of a therapeutic agent or dietary supplement on a cellularfunction of interest. Likewise, models and methods can be used topredict both desirable and undesirable side effects of the therapeuticagent on an interrelated cellular function in the target cell, as wellas the desirable and undesirable effects that may occur in other celltypes. Thus, models and methods disclosed herein can make the drugdevelopment process more rapid and cost effective than is currentlypossible.

The kinetic models can be used to develop individualized drug therapiesand doses. Measuring a set of fluxes and concentrations in an individualbefore and after giving a dose of a therapeutic agent to the individualand constructing kinetic network models according to the measured dataspecific to the individual can be constructed. Analysis of the dynamicand functional properties of the individual specific network can be usedto determine whether the therapeutic agent is effective in the desiredoutcome and the appropriate dose ranges if applicable.

As used herein, the term reaction is intended to mean a conversion thatconsumes a substrate or forms a product that occurs in a biologicalnetwork. The term can include a conversion that occurs due to theactivity of one or more enzymes that are genetically encoded by thehuman genome. The term can also include a conversion that occursspontaneously in a cell, such as in human or mammal. Conversionsincluded in the term include, for example, changes in chemicalcomposition such as those due to nucleophilic or electrophilic addition,nucleophilic or electrophilic substitution, elimination, isomerization,deamination, phosphorylation, methylation, glycolysation, reduction,oxidation or changes in location such as those that occur due to atransport reaction that moves one or more reactants within the samecompartment or from one cellular compartment to another. In the case ofa transport reaction, the substrate and product of the reaction can bechemically the same and the substrate and product can be differentiatedaccording to location in a particular cellular compartment. Thus, areaction that transports a chemically unchanged reactant from a firstcompartment to a second compartment has as its substrate the reactant inthe first compartment and as its product the reactant in the secondcompartment. It will be understood that when used in reference to an insilico model or data structure, a reaction is intended to be arepresentation of a chemical conversion that consumes a substrate orproduces a product.

As used herein, the term reaction is intended to mean a chemical that isa substrate or a product of a reaction that occurs in a biologicalnetwork. The term can include substrates or products of reactionsperformed by one or more enzymes encoded by human gene(s), reactionsoccurring in cells that are performed by one or more non-geneticallyencoded macromolecule, protein or enzyme, or reactions that occurspontaneously in a cell. Metabolites are understood to be reactantswithin the meaning of the term. It will be understood that when used inreference to an in silico model or data structure, a reactant isintended to be a representation of a chemical that is a substrate or aproduct of a reaction that occurs in a cell.

As used herein the term “substrate” is intended to mean a reactant thatcan be converted to one or more products by a reaction. The term caninclude, for example, a reactant that is to be chemically changed due tonucleophilic or electrophilic addition, nucleophilic or electrophilicsubstitution, elimination, isomerization, deamination, phosphorylation,methylation, reduction, oxidation or that is to change location such asby being transported across a membrane or to a different compartment.

As used herein, the term “product” is intended to mean a reactant thatresults from a reaction with one or more substrates. The term caninclude, for example, a reactant that has been chemically changed due tonucleophilic or electrophilic addition, nucleophilic or electrophilicsubstitution, elimination, isomerization, deamination, phosphorylation,methylation, reduction or oxidation or that has changed location such asby being transported across a membrane or to a different compartment.

As used herein, the term “stoichiometric coefficient” is intended tomean a numerical constant correlating the number of one or morereactants and the number of one or more products in a chemical reaction.Typically, the numbers are integers as they denote the number ofmolecules of each reactant in an elementally balanced chemical equationthat describes the corresponding conversion. However, in some cases thenumbers can take on non-integer values, for example, when used in alumped reaction or to reflect empirical data.

As used herein, the term “plurality,” when used in reference toreactions or reactants is intended to mean at least 2 reactions orreactants. The term can include any number of reactions or reactants inthe range from 2 to the number of naturally occurring reactants orreactions for a particular cell. Thus, the term can include, forexample, at least 10, 20, 30, 50, 100, 150, 200, 300, 400, 500, 600 ormore reactions or reactants. The number of reactions or reactants can beexpressed as a portion of the total number of naturally occurringreactions for a particular cell such as at least 20%, 30%, 50%, 60%,75%, 90%, 95% or 98% of the total number of naturally occurringreactions that occur in the particular cell.

As used herein, the term “data structure” is intended to mean a physicalor logical relationship among data elements, designed to supportspecific data manipulation functions. The term can include, for example,a list of data elements that can be added combined or otherwisemanipulated such as a list of representations for reactions from whichreactants can be related in a matrix or network. The term can alsoinclude a matrix that correlates data elements from two or more lists ofinformation such as a matrix that correlates reactants to reactions.Information included in the term can represent, for example, a substrateor product of a chemical reaction, a chemical reaction relating one ormore substrates to one or more products, a constraint placed on areaction, a stoichiometric coefficient, a rate constant, a gradientmatrix coefficient, a gamma matrix coefficient, or a Jacobian matrixcoefficient.

As used herein, the term “constraint”, reaction constraint, or fluxconstraint is intended to mean an upper or lower boundary for areaction. A boundary can specify a minimum or maximum flow of mass,electrons or energy through a reaction. A boundary can further specifydirectionality of a reaction. A boundary can be a constant value such aszero, infinity, or a numerical value such as an integer and non-integer.Alternatively, a boundary can be a variable boundary value as set forthbelow.

As used herein, the term “concentration constraint” is intended to meanan upper and lower boundary for the concentration of a metabolite orchemical species in a biological network. The boundary can specify theamount of mass or number of moles, or number of molecules per unitvolume. This may be within a tissue, a sample of biological fluid, acell, a subcompartment of a cell, or another confined quantity. Theboundary is a non-negative real number including zero, integers, andnon-integers. Alternatively, a boundary for a concentration constraintcan be a variable boundary as set forth below, with the stipulation thatnegative values are not permitted.

As used herein, the term “kinetic constraint” is intended to mean anupper and lower boundary for the kinetic constant for a reaction in abiological network. The boundary can specify the time constantcharacteristic of a particular reaction or set of reactions. Thephysical units are [concentration]^(−n+1)*time⁻¹, where n is the kineticorder of the corresponding reaction. For example, a zero order, firstorder, and second order kinetic constant has units[concentration]*time⁻¹, time⁻¹, and concentration⁻¹*time⁻¹ respectively.The boundary is a non-negative real number including zero, integers, andnon-integers. Alternatively, a boundary for a concentration constraintcan be a variable boundary as set forth below, with the stipulation thatnegative values are not permitted.

As used herein, the term “equilibrium constant constraint” is intendedto mean an upper and lower boundary for the equilibrium constant for areaction in a biological network. The boundary can specify theequilibrium constant of a particular reaction or set of reactions. Theboundary is a non-negative real number including zero, integers, andnon-integers. Alternatively, a boundary for a concentration constraintcan be a variable boundary as set forth below, with the stipulation thatnegative values are not permitted.

As defined herein, the term “uniquely defined constraint” is intended tomean a constraint whose lower bound is equal in value to the upperbound.

As used herein, the term “variable,” when used in reference to aconstraint is intended to mean capable of assuming any of a set ofvalues in response to being acted upon by a constraint function. Theterm “function,” when used in the context of a constraint, is intendedto be consistent with the meaning of the term as it is understood in thecomputer and mathematical arts. A function can be binary such thatchanges correspond to a reaction being off or on. Alternatively,continuous functions can be used such that changes in boundary valuescorrespond to increases or decreases in activity. Such increases ordecreases can also be binned or effectively digitized by a functioncapable of converting sets of values to discreet integer values. Afunction included in the term can correlate a boundary value with thepresence, absence or amount of a biochemical reaction networkparticipant such as a reactant, reaction, enzyme or gene. A functionincluded in the term can correlate a boundary value with an outcome ofat least one reaction in a reaction network that includes the reactionthat is constrained by the boundary limit. A function included in theterm can also correlate a boundary value with an environmental conditionsuch as time, pH, temperature or redox potential.

As used herein, the term “activity,” when used in reference to areaction, is intended to mean the amount of product produced by thereaction, the amount of substrate consumed by the reaction or the rateat which a product is produced or a substrate is consumed. The amount ofproduct produced by the reaction, the amount of substrate consumed bythe reaction or the rate at which a product is produced or a substrateis consumed can also be referred to as the flux for the reaction.

As used herein, the term “constraint set” refers to the collection offlux, concentration, and equilibrium constraints for a biologicalnetwork. Constraint sets for kinetic constants may be found by solving asteady state mass conservation relationship comprising a stoichiometricmatrix and flux vector.

As used herein, the term “allosteric regulation” refers to theregulation of a protein or enzyme by a molecule that binds to a siteother than the primary active site, thereby altering the activity of theenzyme and the corresponding reaction that it regulates.

As used herein, the term “activate” or activation refers to an effect acompound has on another compound, serving to alter the constraints in apositive manner, such as increasing the activity of a reaction. Thisincludes but is not limited to allosteric and non-allosteric regulationof enzymes.

As used herein, the term “inhibit” or inhibition refers to an effect acompound has on another compound, serving to alter the constraints in anegative manner, such as decreasing the activity of a reaction. Thisincludes but is not limited to allosteric and non-allosteric regulationof enzymes.

As used herein, the term “growth” refers to the production of a weightedsum of metabolites identified as biomass components.

As used herein, the term “energy production” refers to the production ofmetabolites that store energy in their chemical bonds, particularly highenergy phosphate bonds such as ATP and GTP.

As used herein, the term “redox equivalent” refers to a metabolite thatcan alter the oxidation state of other metabolites, generally throughthe transfer of electrons and or protons.

As used herein, the term “biomass” refers to a set of metabolites thatare required for cellular functions including replication, maintenance,and/or survival. As such, the term “biomass precursor” refers to ametabolite that can be converted to or is a member of the class ofbiomass metabolites.

As used herein, the term “bioactive small molecule” refers to ametabolite that can inhibit or activate another biological compound,including metabolites, proteins, and nucleic acid polymers.

As used herein, the term “cofactor” refers to a metabolite or chemicalcompound that is required for a catalytic activity of an enzyme to becarried out.

As used herein, the term “optimization problem” refers to a problemwhose solution is obtained through the calculation of a minimum (ormaximum) value of the objective function.

As used herein, the term “objective function” refers to a physiologicalfunction that can be described in terms of the data structure of anetwork, generally in terms of one or more weighted sum or product ofreactions.

Methods and models disclosed herein can be applied to any H. sapienscell type at any stage of differentiation, including, for example,embryonic stem cells, hematopoietic stem cells, differentiatedhematopoietic cells, skeletal muscle cells, cardiac muscle cells, smoothmuscle cells, skin cells, nerve cells, kidney cells, pulmonary cells,liver cells, adipocytes and endocrine cells (e.g., beta islet cells ofthe pancreas, mammary gland cells, adrenal cells, and other specializedhormone secreting cells).

Methods and models disclosed herein can be applied to normal cells orpathological cells. Normal cells that exhibit a variety of physiologicalactivities of interest, including homeostasis, proliferation,differentiation, apoptosis, contraction and motility, can be modeled.Pathological cells can also be modeled, including cells that reflectgenetic or developmental abnormalities, nutritional deficiencies,environmental assaults, infection (such as by bacteria, viral, protozoanor fungal agents), neoplasia, aging, altered immune or endocrinefunction, tissue damage, or any combination of these factors. Thepathological cells can be representative of any type of mammalian,primate or non-primate pathology, including, for example, variousmetabolic disorders of carbohydrate, lipid or protein metabolism,obesity, diabetes, cardiovascular disease, fibrosis, various cancers,kidney failure, immune pathologies, neurodegenerative diseases, andvarious monogenetic diseases described found in humans and documented inthe Online Mendelian Inheritance in Man database (Center for MedicalGenetics, Johns Hopkins University (Baltimore, Md.) and National Centerfor Biotechnology Information, National Library of Medicine (Bethesda,Md.).

Methods and models disclosed herein can also be applied to cellsundergoing therapeutic perturbations, such as cells treated with drugsthat target participants in a reaction network, cells treated withgene-based therapeutics that increase or decrease expression of anencoded protein, and cells treated with radiation. As used herein, theterm “drug” refers to a compound of any molecular nature with a known orproposed therapeutic function, including, for example, small moleculecompounds, peptides and other macromolecules, peptidomimetics andantibodies, any of which can optionally be tagged with cytostatic,targeting or detectable moieties. The term “gene-based therapeutic”refers to nucleic acid therapeutics, including, for example, expressiblegenes with normal or altered protein activity, antisense compounds,ribozymes, DNAzymes, RNA interference compounds (RNAi) and the like. Thetherapeutics can target any reaction network participant, in anycellular location, including participants in extracellular, cellsurface, cytoplasmic, mitochondrial and nuclear locations. Experimentaldata that are gathered on the response of cells to therapeutictreatment, such as alterations in gene or protein expression profiles,can be used to tailor a network for a pathological state of a particularcell type.

Methods and models disclosed herein can be applied to H. sapiens cellsas they exist in any form, such as in primary cell isolates or inestablished cell lines, or in the whole body, in intact organs or intissue explants. Accordingly, the methods and models can take intoaccount intercellular communications and/or inter-organ communications,the effect of adhesion to a substrate or neighboring cells (such as astem cell interacting with mesenchymal cells or a cancer cellinteracting with its tissue microenvironment, or beta-islet cellswithout normal stroma), and other interactions relevant to multicellularsystems.

The reactants to be used in a reaction network data structure of theinvention can be obtained from or stored in a compound database. As usedherein, the term “compound database” is intended to mean a computerreadable medium containing a plurality of molecules that includessubstrates and products of biological reactions. The plurality ofmolecules can include molecules found in multiple organisms, therebyconstituting a universal compound database. Alternatively, the pluralityof molecules can be limited to those that occur in a particularorganism, thereby constituting an organism-specific compound database.Each reactant in a compound database can be identified according to thechemical species and the cellular compartment in which it is present.Thus, for example, a distinction can be made between glucose in theextracellular compartment versus glucose in the cytosol. Additionallyeach of the reactants can be specified as a metabolite of a primary orsecondary metabolic pathway. Although identification of a reactant as ametabolite of a primary or secondary metabolic pathway does not indicateany chemical distinction between the reactants in a reaction, such adesignation can assist in visual representations of large networks ofreactions.

As used herein, the term “compartment” is intended to mean a subdividedregion containing at least one reactant, such that the reactant isseparated from at least one other reactant in a second region. Asubdivided region included in the term can be correlated with asubdivided region of a cell. Thus, a subdivided region included in theterm can be, for example, the intracellular space of a cell; theextracellular space around a cell; the periplasmic space; the interiorspace of an organelle such as a mitochondrium, endoplasmic reticulum,Golgi apparatus, vacuole or nucleus; or any subcellular space that isseparated from another by a membrane or other physical barrier.Subdivided regions can also be made in order to create virtualboundaries in a reaction network that are not correlated with physicalbarriers. Virtual boundaries can be made for the purpose of segmentingthe reactions in a network into different compartments or substructures.

As used herein, the term “substructure” is intended to mean a portion ofthe information in a data structure that is separated from otherinformation in the data structure such that the portion of informationcan be separately manipulated or analyzed. The term can include portionssubdivided according to a biological function including, for example,information relevant to a particular metabolic pathway such as aninternal flux pathway, exchange flux pathway, central metabolic pathway,peripheral metabolic pathway, or secondary metabolic pathway. The termcan include portions subdivided according to computational ormathematical principles that allow for a particular type of analysis ormanipulation of the data structure.

The reactions included in a reaction network data structure can beobtained from a metabolic reaction database that includes thesubstrates, products, and stoichiometry of a plurality of biologicalreactions. The reactants in a reaction network data structure can bedesignated as either substrates or products of a particular reaction,each with a stoichiometric coefficient assigned to it to describe thechemical conversion taking place in the reaction. Each reaction is alsodescribed as occurring in either a reversible or irreversible direction.Reversible reactions can either be represented as one reaction thatoperates in both the forward and reverse direction or be decomposed intotwo irreversible reactions, one corresponding to the forward reactionand the other corresponding to the backward reaction.

Reactions included in a reaction network data structure can includeintra-system or exchange reactions. Intra-system reactions are thechemically and electrically balanced interconversions of chemicalspecies and transport processes, which serve to replenish or drain therelative amounts of certain metabolites. These intra-system reactionscan be classified as either being transformations or translocations. Atransformation is a reaction that contains distinct sets of compounds assubstrates and products, while a translocation contains reactantslocated in different compartments. Thus, a reaction that simplytransports a metabolite from the extracellular environment to thecytosol, without changing its chemical composition is solely classifiedas a translocation, while a reaction such as the phosphotransferasesystem (PTS) which takes extracellular glucose and converts it intocytosolic glucose-6-phosphate is a translocation and a transformation.

Exchange reactions are those which constitute sources and sinks,allowing the passage of metabolites into and out of a compartment oracross a hypothetical system boundary. These reactions are included in amodel for simulation purposes and represent the metabolic demands placedon a cell. While they may be chemically balanced in certain cases, theyare typically not balanced and can often have only a single substrate orproduct. As a matter of convention the exchange reactions are furtherclassified into demand exchange and input/output exchange reactions.

Input/output exchange reactions are used to allow extracellularreactants to enter or exit the reaction network represented by a modelof the invention. For each of the extracellular metabolites acorresponding input/output exchange reaction can be created. Thesereactions can either be irreversible or reversible with the metaboliteindicated as a substrate with a stoichiometric coefficient of one and noproducts produced by the reaction. This particular convention is adoptedto allow the reaction to take on a positive flux value (activity level)when the metabolite is being produced or removed from the reactionnetwork and a negative flux value when the metabolite is being consumedor introduced into the reaction network. These reactions will be furtherconstrained during the course of a simulation to specify exactly whichmetabolites are available to the cell and which can be excreted by thecell.

A demand exchange reaction is always specified as an irreversiblereaction containing at least one substrate. These reactions aretypically formulated to represent the production of an intracellularmetabolite by the metabolic network or the aggregate production of manyreactants in balanced ratios such as in the representation of a reactionthat leads to biomass formation, also referred to as growth.

A demand exchange reaction can be introduced for any metabolite in amodel of the invention. Most commonly these reactions are introduced formetabolites that are required to be produced by the cell for thepurposes of creating a new cell such as amino acids, nucleotides,phospholipids, and other biomass constituents, or metabolites that areto be produced for alternative purposes. Once these metabolites areidentified, a demand exchange reaction that is irreversible andspecifies the metabolite as a substrate with a stoichiometriccoefficient of unity can be created. With these specifications, if thereaction is active it leads to the net production of the metabolite bythe system meeting potential production demands. Examples of processesthat can be represented as a demand exchange reaction in a reactionnetwork data structure and analyzed by the methods of the inventioninclude, for example, production or secretion of an individual protein;production or secretion of an individual metabolite such as an aminoacid, vitamin, nucleoside, antibiotic or surfactant; production of ATPfor extraneous energy requiring processes such as locomotion; orformation of biomass constituents.

In addition to these demand exchange reactions that are placed onindividual metabolites, demand exchange reactions that utilize multiplemetabolites in defined stoichiometric ratios can be introduced. Thesereactions are referred to as aggregate demand exchange reactions. Anexample of an aggregate demand reaction is a reaction used to simulatethe concurrent growth demands or production requirements associated withcell growth that are placed on a cell, for example, by simulating theformation of multiple biomass constituents simultaneously at aparticular cellular growth rate.

Depending upon the particular environmental conditions being tested andthe desired activity, a reaction network data structure can containsmaller numbers of reactions such as at least 200, 150, 100 or 50reactions. A reaction network data structure having relatively fewreactions can provide the advantage of reducing computation time andresources required to perform a simulation. When desired, a reactionnetwork data structure having a particular subset of reactions can bemade or used in which reactions that are not relevant to the particularsimulation are omitted. Alternatively, larger numbers of reactions canbe included in order to increase the accuracy or molecular detail of themethods of the invention or to suit a particular application. Thus, areaction network data structure can contain at least 300, 350, 400, 450,500, 550, 600 or more reactions up to the number of reactions that occurin a particular cell or that are desired to simulate the activity of thefull set of reactions occurring in the particular human cell.

A reaction network data structure or index of reactions used in the datastructure such as that available in a metabolic reaction database, asdescribed herein, can be annotated to include information about aparticular reaction. A reaction can be annotated to indicate, forexample, assignment of the reaction to a protein, macromolecule orenzyme that performs the reaction, assignment of a gene(s) that codesfor the protein, macromolecule or enzyme, the Enzyme Commission (EC)number of the particular metabolic reaction or Gene Ontology (GO) numberof the particular metabolic reaction, a subset of reactions to which thereaction belongs, citations to references from which information wasobtained, or a level of confidence with which a reaction is believed tooccur in a particular human cell. A computer readable medium of theinvention can include a gene database containing annotated reactions.Such information can be obtained during the course of building ametabolic reaction database or model of the invention as describedbelow.

Thus, a method is provided for making or expanding a data structurerelating a plurality of H. sapiens reactants to a plurality of H.sapiens reactions in a computer readable medium. The method includes thesteps of: (a) identifying a plurality of H. sapiens reactions and aplurality of H. sapiens reactants that are substrates and products ofthe H. sapiens reactions; (b) relating the plurality of H. sapiensreactants to the plurality of H. sapiens reactions in a data structure,wherein each of the H. sapiens reactions includes one or more reactantsidentified as a substrate of the reaction, one or more reactantsidentified as a product of the reaction and a stoichiometric coefficientrelating the substrate and the product; (c) making a constraint set forthe plurality of H. sapiens reactions; (d) providing an objectivefunction; (e) determining at least one flux distribution that minimizesor maximizes the objective function when the constraint set is appliedto the data structure, (f) if at least one flux distribution is notpredictive of human physiology, then adding a reaction to or deleting areaction from the data structure and repeating step (e), and (g) if atleast one flux distribution is predictive of human physiology, thenstoring the data structure in a computer readable medium.

The majority of the reactions occurring in human reaction networks arecatalyzed by enzymes/proteins, which are created through thetranscription and translation of the genes found on the chromosome(s) inthe cell. The remaining reactions occur through non-enzymatic processes.Furthermore, a reaction network data structure can contain reactionsthat add or delete steps to or from a particular reaction pathway. Forexample, reactions can be added to optimize or improve performance of aH. sapiens model in view of empirically observed activity.Alternatively, reactions can be deleted to remove intermediate steps ina pathway when the intermediate steps are not necessary to model fluxthrough the pathway. For example, if a pathway contains threenonbranched steps, the reactions can be combined or added together togive a net reaction, thereby reducing memory required to store thereaction network data structure and the computational resources requiredfor manipulation of the data structure.

A reaction network data structure as described herein can be used todetermine the activity of one or more reactions in a plurality of H.sapiens reactions independent of any knowledge or annotation of theidentity of the protein that performs the reaction or the gene encodingthe protein. A model that is annotated with gene or protein identitiescan include reactions for which a protein or encoding gene is notassigned. While a large portion of the reactions in a cellular metabolicnetwork are associated with genes in the organism's genome, there arealso a substantial number of reactions included in a model for whichthere are no known genetic associations. Such reactions can be added toa reaction database based upon other information that is not necessarilyrelated to genetics such as biochemical or cell based measurements ortheoretical considerations based on observed biochemical or cellularactivity. For example, there are many reactions that are notprotein-enabled reactions. Furthermore, the occurrence of a particularreaction in a cell for which no associated proteins or genetics havebeen currently identified can be indicated during the course of modelbuilding by the iterative model building methods as disclosed herein.

The reactions in a reaction network data structure or reaction databasecan be assigned to subsystems by annotation, if desired. The reactionscan be subdivided according to biological criteria, such as according totraditionally identified metabolic pathways (glycolysis, amino acidmetabolism and the like) or according to mathematical or computationalcriteria that facilitate manipulation of a model that incorporates ormanipulates the reactions. Methods and criteria for subdividing areaction database are described in further detail in Schilling et al.,J. Theor. Biol. 203:249-283 (2000), which is hereby incorporated byreference in its entirety. The use of subsystems can be advantageous fora number of analysis methods, such as extreme pathway analysis, and canmake the management of model content easier. Although assigningreactions to subsystems can be achieved without affecting the use of theentire model for simulation, assigning reactions to subsystems can allowa user to search for reactions in a particular subsystem, which may beuseful in performing various types of analyses. Therefore, a reactionnetwork data structure can include any number of desired subsystemsincluding, for example, 2 or more subsystems, 5 or more subsystems, 10or more subsystems, 25 or more subsystems or 50 or more subsystems.

Flux constraints can be placed on the value of any of the fluxes in themetabolic network using a constraint set. These constraints can berepresentative of a minimum or maximum allowable flux through a givenreaction, possibly resulting from a limited amount of an enzyme present.Additionally, the constraints can determine the direction orreversibility of any of the reactions or transport fluxes in thereaction network data structure.

The ability of a reaction to be actively occurring is dependent on alarge number of additional factors beyond just the availability ofsubstrates. These factors, which can be represented as variableconstraints in the models and methods of the invention include, forexample, the presence of cofactors necessary to stabilize theprotein/enzyme, the presence or absence of enzymatic inhibition andactivation factors, the active formation of the protein/enzyme throughtranslation of the corresponding mRNA transcript, the transcription ofthe associated gene(s) or the presence of chemical signals and/orproteins that assist in controlling these processes that ultimatelydetermine whether a chemical reaction is capable of being carried outwithin an organism. Of particular importance in the regulation of humancell types is the implementation of paracrine and endocrine signalingpathways to control cellular activities. In these cases a cell secretessignaling molecules that may be carried far afield to act on distanttargets (endocrine signaling), or act as local mediators (paracrinesignaling). Examples of endocrine signaling molecules include hormonessuch as insulin, while examples of paracrine signaling molecules includeneurotransmitters such as acetylcholine. These molecules induce cellularresponses through signaling cascades that affect the activity ofbiochemical reactions in the cell. Regulation can be represented in anin silico H. sapiens model by providing a variable constraint.

The methods described herein can be implemented on any conventional hostcomputer system, such as those based on Intel® or AMID® microprocessorsand running Microsoft Windows operating systems. Other systems, such asthose using the UNIX or LINUX operating system and based on IBM®, DEC®or Motorola® microprocessors are also contemplated. The systems andmethods described herein can also be implemented to run on client-serversystems and wide-area networks, such as the Internet.

Software to implement a method or model of the invention can be writtenin any well-known computer language, such as Java, C, C++, Visual Basic,FORTRAN or COBOL and compiled using any well-known compatible compiler.The software of the invention normally runs from instructions stored ina memory on a host computer system. A memory or computer readable mediumcan be a hard disk, floppy disc, compact disc, DVD, magneto-opticaldisc, Random Access Memory, Read Only Memory or Flash Memory. The memoryor computer readable medium used in the invention can be containedwithin a single computer or distributed in a network. A network can beany of a number of conventional network systems known in the art such asa local area network (LAN) or a wide area network (WAN). Client-serverenvironments, database servers and networks that can be used in theinvention are well known in the art. For example, the database servercan run on an operating system such as UNIX, running a relationaldatabase management system, a World Wide Web application and a WorldWide Web server. Other types of memories and computer readable media arealso contemplated to function within the scope of the invention.

The data matrices constructed by the methods described in this inventioncan be represented in a markup language format including, for example,Standard Generalized Markup Language (SGML), Hypertext markup language(HTML) or Extensible Markup language (XML). Markup languages can be usedto tag the information stored in a database or data structure of theinvention, thereby providing convenient annotation and transfer of databetween databases and data structures. In particular, an XML format canbe useful for structuring the data representation of reactions,reactants and their annotations; for exchanging database contents, forexample, over a network or internet; for updating individual elementsusing the document object model; or for providing differential access tomultiple users for different information content of a data base or datastructure of the invention. XML programming methods and editors forwriting XML code are known in the art as described, for example, in Ray,Learning XML O'Reilly and Associates, Sebastopol, Calif. (2001).

A computer system of the invention can further include a user interfacecapable of receiving a representation of one or more reactions. A userinterface of the invention can also be capable of sending at least onecommand for modifying the data structure, the constraint set or thecommands for applying the constraint set to the data representation, ora combination thereof. The interface can be a graphic user interfacehaving graphical means for making selections such as menus or dialogboxes. The interface can be arranged with layered screens accessible bymaking selections from a main screen. The user interface can provideaccess to other databases useful in the invention such as a metabolicreaction database or links to other databases having informationrelevant to the reactions or reactants in the reaction network datastructure or to human physiology. Also, the user interface can display agraphical representation of a reaction network or the results of asimulation using a model of the invention.

Once an initial reaction network data structure and set of constraintshas been created, a model disclosed herein can be tested by preliminarysimulation. During preliminary simulation, gaps in the network or“dead-ends” in which a metabolite can be produced but not consumed orwhere a metabolite can be consumed but not produced can be identified.Based on the results of preliminary simulations areas of the metabolicreconstruction that require an additional reaction can be identified.The determination of these gaps can be readily calculated throughappropriate queries of the reaction network data structure and need notrequire the use of simulation strategies, however, simulation would bean alternative approach to locating such gaps.

In the preliminary simulation testing and model content refinement stagean existing model may be subjected to a series of functional tests todetermine if it can perform basic requirements such as the ability toproduce the required biomass constituents and generate predictionsconcerning the basic physiological characteristics of the particularorganism strain being modeled. The more preliminary testing that isconducted, typically, the higher the quality of the model that will begenerated. Typically the majority of the simulations used in this stageof development will be single optimizations. A single optimization canbe used to calculate a single flux distribution demonstrating howmetabolic resources are routed determined from the solution to oneoptimization problem. An optimization problem can be solved using linearprogramming as demonstrated in the Examples below. The result can beviewed as a display of a flux distribution on a reaction map. Temporaryreactions can be added to the network to determine if they should beincluded into the model based on modeling/simulation requirements.

Once a model is sufficiently complete with respect to the content of thereaction network data structure according to the criteria set forthabove, the model can be used to simulate activity of one or morereactions in a reaction network. The results of a simulation can bedisplayed in a variety of formats including, for example, a table,graph, reaction network, flux distribution map or a as a modal matrix.

FIG. 6 is a block diagram of some contemplated therapeutic anddiagnostic embodiments of the invention employing the model. Diagnosticor therapeutic analysis 600 begins 601 by first gathering the archetypedata 602 for the domain under investigation and constructing S usingthis data. One skilled in the art will readily recognize numerous datasets that are available for a given domain. For example, in the domainof drug therapy identification, archetypical data may comprisestatistical and theoretical information regarding the reactions underinvestigation. Physical and context specific data is then incorporated603. Physical and context specific data may comprise the metabolomiccharacteristics common to a patient demographic under investigation. Themodel is then used to simulate steady state conditions 604. Once thesteady state is achieved, metabolic data and equilibrium constant dataare incorporated 605. The steady state equations are again solved for,and the matrix G and the dynamic rate equations are derived 606. At thispoint, the dynamic model provides an analysis of Jx and Jv 607 orsimulations of the rate equations 608. The modeled dynamics are observed609 and features of interest identified 610. One skilled in the art willrecognize that features of interest will depend on the domain underinvestigation and metabolic pathway being studied. For example, fordiagnostic purposes, particular fluxes or concentrations whose variationproduces novel and unexpected behavior in the metabolic pathway may beconsidered features of interest.

If no features are observed, or further analysis is required, parametersor steady state conditions are modified in the model 611 so as tosimulate a different dynamic, and the simulation 604 observed againuntil features of interest may be identified.

Alternatively, where features of interest have been identified,physiological materials relevant to them may be prepared 612 and theprocess completes 613. In the domain of drug therapy, for example, thepreparation of physiological materials may comprise the identificationand synthesis of compounds which will generate the desired metabolicbehavior responsible for generating the features of interest previouslymodeled. One skilled in the art, however, will readily recognize thatthe model will lend itself to a wide variety of domains requiring acomprehensive understanding of large and complex metabolic pathways.Accordingly, the information derived from the model may be used innumerous physical applications, from confirming disease diagnosis, toidentifying environmental factors impacting individual health.

Thus, methods provided herein can predict physiological function. In oneembodiment, a method is provided, the method including providing a datastructure relating a plurality of H. sapiens reactants to a plurality ofH. sapiens reactions, wherein each of the H. sapiens reactions includesone or more reactants identified as a substrate of the reaction, one ormore reactants identified as a product of the reaction and astoichiometric coefficient relating the substrate and the product;providing a constraint set for the plurality of H. sapiens reactions;providing an objective function, and determining at least one fluxdistribution that minimizes or maximizes the objective function when theconstraint set is applied to the data structure, thereby predicting a H.sapiens physiological function.

As used herein, the term “physiological function,” when used inreference to H. sapiens, is intended to mean an activity of a human cellas a whole. An activity included in the term can be the magnitude orrate of a change from an initial state of a human cell to a final stateof the human cell. An activity can be measured qualitatively orquantitatively. An activity included in the term can be, for example,growth, energy production, redox equivalent production, biomassproduction, development, or consumption of carbon, nitrogen, sulfur,phosphate, hydrogen or oxygen. An activity can also be an output of aparticular reaction that is determined or predicted in the context ofsubstantially all of the reactions that affect the particular reactionin a human cell or substantially all of the reactions that occur in ahuman cell. Examples of a particular reaction included in the term areproduction of biomass precursors, production of a protein, production ofan amino acid, production of a purine, production of a pyrimidine,production of a lipid, production of a fatty acid, production of acofactor, production of a hormone, production of a bioactive smallmolecule, or transport of a metabolite. A physiological function caninclude an emergent property which emerges from the whole but not fromthe sum of parts where the parts are observed in isolation (see forexample, Palsson Nat. Biotech 18:1147-1150 (2000)).

A physiological function of H. sapiens reactions can be determined usingdynamic phase plane analysis.

A physiological function of H. sapiens reactions can also be determinedusing a reaction map to display a flux distribution. A reaction map canbe used to view reaction networks at a variety of levels. In the case ofa cellular metabolic reaction network, a reaction map can contain theentire reaction complement representing a global perspective.Alternatively, a reaction map can focus on a particular region ofmetabolism such as a region corresponding to a reaction subsystemdescribed above or even on an individual pathway or reaction.

Methods disclosed herein can be used to determine the activity of aplurality of H. sapiens reactions including, for example, biosynthesisof an amino acid, degradation of an amino acid, biosynthesis of apurine, biosynthesis of a pyrimidine, biosynthesis of a lipid,metabolism of a fatty acid, biosynthesis of a cofactor, production of ahormone, production of a bioactive small molecule, transport of ametabolite and metabolism of an alternative carbon source.

Methods disclosed herein can be used to determine a phenotype of a H.sapiens mutant. The activity of one or more H. sapiens reactions can bedetermined using the methods described above, wherein the reactionnetwork data structure lacks one or more gene-associated reactions thatoccur in H. sapiens. Alternatively, methods can be used to determine theactivity of one or more H. sapiens reactions when a reaction that doesnot naturally occur in H. sapiens is added to the reaction network datastructure. Deletion of a gene can also be represented in a model of theinvention by constraining the flux through the reaction to zero, therebyallowing the reaction to remain within the data structure. Thus,simulations can be made to predict the effects of adding or removinggenes to or from H. sapiens. The methods can be particularly useful fordetermining the effects of adding or deleting a gene that encodes for agene product that performs a reaction in a peripheral metabolic pathway.In one contemplated embodiment, adenovirus vectors are used for in vivotransfer of genes determined in silico to be required for a desiredfunctioning of the metabolic pathway.

A drug target or target for any other agent that affects human cellularfunction can be predicted using the methods of the invention. Suchpredictions can be made by removing a reaction to simulate totalinhibition or prevention by a drug or agent. Alternatively, partialinhibition or reduction in the activity a particular reaction can bepredicted by performing the methods with altered constraints. Themethods can be particularly useful for identifying a target in aperipheral metabolic pathway.

Once a reaction has been identified for which activation or inhibitionproduces a desired effect on H. sapiens function, an enzyme ormacromolecule that performs the reaction in H. sapiens or a gene thatexpresses the enzyme or macromolecule can be identified as a target fora drug or other agent. A candidate compound for a target identified bythe methods of the invention can be isolated or synthesized using knownmethods. Such methods for isolating or synthesizing compounds caninclude, for example, rational design based on known properties of thetarget (see, for example, Decamp et al., Protein Engineering Principlesand Practice, Ed. Cleland and Craik, Wiley-Liss, New York, pp. 467-506(1996), which is hereby incorporated by reference in its entirety),screening the target against combinatorial libraries of compounds (seefor example, Houghten et al., Nature, 354, 84-'86 (1991); Dooley et al.,Science, 266, 2019-2022 (1994), which describe an iterative approach, orR. Houghten et al. PCT/US91/08694 and U.S. Pat. No. 5,556,762 whichdescribe the positional scanning approach, all four of which are herebyincorporated by reference in their entireties), or a combination of bothto obtain focused libraries. Those skilled in the art will know or willbe able to routinely determine assay conditions to be used in a screenbased on properties of the target or activity assays known in the art.

A candidate drug or agent, whether identified by methods described aboveor by other methods known in the art, can be validated using an insilico human metabolic model or method disclosed herein. The effect of acandidate drug or agent on human physiological function can be predictedbased on the activity for a target in the presence of the candidate drugor agent measured in vitro or in vivo. This activity can be representedin an in silico human metabolic model by adding a reaction to the model,removing a reaction from the model or adjusting a constraint for areaction in the model to reflect the measured effect of the candidatedrug or agent on the activity of the reaction. By running a simulationunder these conditions the holistic effect of the candidate drug oragent on human physiological function can be predicted.

Once a candidate has been identified, medical practitioners may readilyconceive a plurality of delivery methods that will induce the desiredsimulated metabolic change. Suitable pharmaceutically acceptablecarriers and their formulation are described in standard formulationtreatises, e.g. Remington's Pharmaceuticals Sciences by E. W. martin.See also Wang, Y. J. and Hanson, M. A. “Parental formulations ofProteins and Pepties; Stability and Stabilizers,” Journals of ParentalSciences and Technology, Technical Report No. 10, Supp. 42:2 S (1988).

The methods of the invention can be used to determine the effects of oneor more environmental components or conditions on an activity of a H.sapiens cell. As set forth above, an exchange reaction can be added to areaction network data structure corresponding to uptake of anenvironmental component, release of a component to the environment, orother environmental demand. The effect of the environmental component orcondition can be further investigated by running simulations withadjusted values for the metabolic flux vector of the exchange reactiontarget reaction to reflect a finite maximum or minimum flux valuecorresponding to the effect of the environmental component or condition.The environmental component can be, for example an alternative carbonsource or a metabolite that when added to the environment of a H.sapiens cell can be taken up and metabolized. The environmentalcomponent can also be a combination of components present for example ina minimal medium composition. Thus, methods disclosed herein can be usedto determine an optimal or minimal medium composition that is capable ofsupporting a particular activity of H. sapiens.

In some embodiments, a method for determining a set of environmentalcomponents to achieve a desired activity for H. sapiens is provided. Themethod includes the steps of providing a data structure relating aplurality of H. sapiens reactants to a plurality of H. sapiensreactions, wherein each of the H. sapiens reactions includes one or morereactants identified as a substrate of the reaction, one or morereactants identified as a product of the reaction and a stoichiometriccoefficient relating the substrate and the product; providing aconstraint set for the plurality of H. sapiens reactions; applying theconstraint set to the data representation, thereby determining theactivity of one or more H. sapiens reactions; determining the activityof one or more Homo sapiens reactions according to steps (a) through(c), wherein the constraint set includes an upper or lower bound on theamount of an environmental component, and repeating steps (a) through(c) with a changed constraint set, wherein the activity determined.

Further contemplated embodiments comprise down-regulating portions ofmetabolic networks by introducing inhibiting compounds regulating flowsfound in the model. These compounds may inhibit deleterious feedbacksystems, such as by suppressing antigen processing or antigenpresentation. Other embodiments may up-regulate portions of metabolicnetworks to account for deletions or omissions in preceding stages of a“normally” functioning network. Such pathways need not be confined to H.sapiens but can be used in foreign entities as well. For example,metabolic pathways in viruses and bacteria can be modeled to determinean exogenous chemical dependency which will permit the virus' liveintroduction into a host, and subsequent controlled attenuation viareduction of the chemical dependency. Such live introduction facilitatesbetter immune response. Embodiments of the present invention, bymodeling the metabolic pathway, permit more controlled attenuation ofthese vectors, with less risk to the host.

One skilled in the art will readily recognize that these are but a fewpossible implementations of the invention and that many more exist.

The following examples are intended to illustrate but not limit theinvention.

Example 1 A Stoichiometric Model of Glycolysis in Homo sapiens

The stoichiometric matrix as described in this invention was constructedfor the glycolytic pathway in Homo sapiens.

TABLE 1B The stoichiometric matrix for the glycolytic pathway, thecolumns correspond to metabolites and columns correspond to fluxes. HKPGI PFK ALD TPI GAPDH PGK PGM EN PK LDH LEX ATPase G6P 1 −1 0 0 0 0 0 00 0 0 0 0 F6P 0 1 −1 0 0 0 0 0 0 0 0 0 0 FDP 0 0 1 −1 0 0 0 0 0 0 0 0 0DHAP 0 0 0 1 −1 0 0 0 0 0 0 0 0 GAP 0 0 0 1 1 −1 0 0 0 0 0 0 0 DPG13 0 00 0 0 1 −1 0 0 0 0 0 0 PG3 0 0 0 0 0 0 1 −1 0 0 0 0 0 PG2 0 0 0 0 0 0 01 −1 0 0 0 0 PEP 0 0 0 0 0 0 0 0 1 −1 0 0 0 PYR 0 0 0 0 0 0 0 0 0 1 −1 00 LAC 0 0 0 0 0 0 0 0 0 0 1 −1 0 ATP −1 0 −1 0 0 0 1 0 0 1 0 0 −1

Example 2 The Gradient Matrix of Glycolysis in Homo sapiens

A model of glycolysis was constructed by the methods described in thisinvention. Fluxes and concentrations for a normal human red blood cellare listed in Table 2. Table 2 also includes equilibrium constants forthe reactions in the network. Reactions with equilibrium constants muchgreater than 1 or much less than one are effectively irreversible. Thegradient matrix can be factored into the kappa and gamma matrices asdescribed in this invention and illustrated for the glycolytic pathwayin Homo sapiens.

TABLE 2 The steady state concentrations of metabolites (mM: millimolar).Concentrations (mM) G6P 0.066440083 F6P 0.027233463 FDP 0.013861036 DHAP0.147408623 GAP 0.006698582 DPG13 0.000228927 PG3 0.083927674 PG20.012267539 PEP 0.020874113 PYR 0.06 LAC 1.397092049 ATP 2.401613356

TABLE 3 The steady state fluxes and the equilibrium constants for thereactions in the glycolytic network (mM: millimolar, hr: hour, Keq:equilibrium constant). Flux (mM/hr) Keq HK 0.498225636 1000 PGI0.498225611 0.41 PFK 0.498225803 1000 ALD 0.498225611 12.34568 TPI0.498225611 0.057143 GAPDH 0.996253719 55.86592 PGK −0.714792222 1800PGM 0.996451222 0.147059 EN 0.996451222 1000 PK 0.99270572 1000 LDH13.71530938 4464.286 LEX 0.996451222 1 ATPase 0.996451975 1000

By solving the steady state mass balance equations using thestoichiometric matrix in Example 1 and the mass-action rate lawformulation and substituting in steady state metabolite concentrations(Table 2), steady state flux concentrations (Table 3), and equilibriumconstants (Table 3), the forward rate constants were calculated.

The mass-action rate laws were then written in terms of the equilibriumconstants and rate constants, as function of the metaboliteconcentrations. The gradient matrix was calculated from this.

TABLE 4 The gradient matrix for the glycolytic pathway. Each element inthe row is the partial derivative of the corresponding reaction withrespect to the metabolite in the corresponding column. G6P F6P FDP DHAPGAP DPG13 PG3 HK   −0.000101     0  0  0    0    0     0 PGI 29304.9−71475.36  0  0    0    0     0 PFK   0     18.29651 −0.003721  0    0   0     0 ALD   0     0 36.15294 −0.019616  −0.43167    0     0 TPI   0    0  0 16.50659 −288.8654    0     0 GAPDH   0     0  0  0   148.7583   −0.94482   −0.00061 PGK   0     0  0  0    0 2009588 −5490.018 PGM  0     0  0  0    0    0   1959.931 EN   0     0  0  0    0    0     0PK   0     0  0  0    0    0     0 LDH   0     0  0  0    0    39.25606   39.25606 LEX   0     0  0  0    0    0     0 ATPase   0     0  0  0   0    0     0 PG2 PEP PYR LAC ATP HK     0  0  0   0     0.207471 PGI    0  0  0   0     0 PFK     0  0  0   0     0.207582 ALD     0  0  0  0     0 TPI     0  0  0   0     0 GAPDH   −0.00061 −0.00061  −0.00061  0     0 PGK     0  0  0   0 −1133.834 PGM −13327.53  0  0   0     0 EN    81.36511 −0.081365  0   0     0 PK     0 48.23862  −0.237211   0  −2.067691 LDH     39.25606 39.25606 271.2571 −0.146556     0 LEX     0 0  0   5.055766     0 ATPase     0  0  0   0     0.415409

The structural similarity between the stoichiometric matrix and negativeof the transpose of the gradient matrix for glycolysis is immediatelyapparent by comparing Table 1 and Table 4.

As discussed above, the gradient matrix can be factored into the kappaand gamma matrices as illustrated in Table 5 and Table 6.

TABLE 5 The kappa matrix which was constructed by factoring out thenorms in the gradient matrix in Table 4. The 1-norm was used for thisglycolytic kappa matrix. HK 0.207572    0  0  0  0  0    0 PGI 0100780.3  0  0  0  0    0 PFK 0    0 18.50781  0  0  0    0 ALD 0    0 0 36.60422  0  0    0 TPI 0    0  0  0 305.372  0    0 GAPDH 0    0  0 0  0 149.7056    0 PGK 0    0  0  0  0  0 2016211 PGM 0    0  0  0  0 0    0 EN 0    0  0  0  0  0    0 PK 0    0  0  0  0  0    0 LDH 0    0 0  0  0  0    0 LEX 0    0  0  0  0  0    0 ATPase 0    0  0  0  0  0   0 HK   0  0  0  0 0 0 PGI   0  0  0  0 0 0 PFK   0  0  0  0 0 0 ALD  0  0  0  0 0 0 TPI   0  0  0  0 0 0 GAPDH   0  0  0  0 0 0 PGK   0  0 0  0 0 0 PGM 15287.46  0  0  0 0 0 EN   0 81.44648  0  0 0 0 PK   0  050.54352  0 0 0 LDH   0  0  0 428.4279 0 0 LEX   0  0  0  0 5.055766 0ATPase   0  0  0  0 0 0.415409

TABLE 6 The gamma matrix corresponding to the kappa matrix in Table 5and the gradient matrix in Table 4. G6P F6P FDP DHAP GAP DPG13 PG3−0.000488   0   0   0   0   0   0   0.29078 −0.70922   0   0   0   0   0  0   0.988583 −0.000201   0   0   0   0   0   0   0.987671 −0.000536−0.011793   0   0   0   0   0   0.054054 −0.945946   0   0   0   0   0  0   0.993673 −0.006311 −4.07E−06   0   0   0   0   0   0.996715−0.002723   0   0   0   0   0   0   0.128205   0   0   0   0   0   0   0  0   0   0   0   0   0   0   0   0   0   0   0   0.091628   0.091628  0   0   0   0   0   0   0   0   0   0   0   0   0   0 PG2 PEP PYR LACATP   0   0   0   0   0.999512   0   0   0   0   0   0   0   0   0  0.011216   0   0   0   0   0   0   0   0   0   0 −4.07E−06 −4.07E−06−4.07E−06   0   0   0   0   0   0 −0.000562 −0.871795   0   0   0   0  0.999001 −0.000999   0   0   0   0   0.954398 −0.004693   0 −0.040909  0.091628   0.091628   0.633145 −0.000342   0   0   0   0   1   0   0  0   0   0   1

The elements in the kappa matrix are spread over approximately ten (log[κ_(max)/κ_(min)]=9.7) orders of magnitude, which conveys thebiologically important time scale hierarchy and also conveys theill-conditioned properties of the gradient matrix.

Example 3 The Dynamic Glycolytic Pathway in Homo sapiens

Using the stoichiometric matrix in Example 1 and the gradient matrix inExample 2, a dynamic model of the glycolytic pathway was constructed bycalculating the concentration Jacobian matrix by post-multiplication ofthe stoichiometric matrix by the gradient matrix (Table 7) and the fluxJacobian matrix by pre-multiplication of the stoichiometric matrix bythe gradient matrix (Table 8), as per the methods described above.

TABLE 7 The concentration Jacobian matrix for the glycolytic pathway inHomo sapiens. Row and column headings are labeled for clarity. G6P F6PFDP DHAP GAP DPG13 PG3 G6P −29304.9   71475.36    0    0    0      0    0 F6P   29304.9 −71493.65    0.003721    0    0      0     0 FDP    0     18.29651 −36.15666    0.019616    0.43167      0     0 DHAP    0     0   36.15294 −16.52621   288.4337      0     0 GAP     0     0  36.15294   16.48698 −438.0554      0.94482     0.00061 DPG13     0    0    0    0   148.7583 −2009588   5490.017 PG3     0     0    0    0   0   2009588 −7449.948 PG2     0     0    0    0    0      0  1959.931 PEP     0     0    0    0    0      0     0 PYR     0     0   0    0    0    −39.25606  −39.25606 LAC     0     0    0    0    0     39.25606    39.25606 ATP     0.000101   −18.29651    0.003721    0   0   2009588 −5490.018 PG2 PEP PYR LAC ATP G6P     0    0    0   0    0.207471 F6P     0    0    0   0   −0.207582 FDP     0    0    0   0    0.207582 DHAP     0    0    0   0     0 GAP     0.00061    0.00061   0.00061   0     0 DPG13   −0.00061  −0.00061  −0.00061   0   1133.834PG3   13327.53    0    0   0 −1133.834 PG2 −13408.89    0.081365    0  0     0 PEP     81.36511 −48.31998    0.237211   0     2.067691 PYR  −39.25606    8.982558 −271.4943   0.146556   −2.067691 LAC    39.25606   39.25606   271.2571 −5.202322     0 ATP     0   48.23862 −0.237211   0 −1136.733

TABLE 8 The concentration Jacobian matrix for the glycolytic pathway inHomo sapiens. Row and column headings are labeled for clarity. HK PGIPFK ALD TPI GAPDH PGK HK   −0.207572      0.000101   −0.207471    0    0   0      0.207471 PGI 29304.9 −100780.3 71475.36    0    0    0      0PFK   −0.207582     18.29651  −18.50781    0.003721    0    0     0.207582 ALD   0      0   36.15294  −36.60422  −0.412054    0     0 TPI   0      0   0 −272.3588 −305.372   288.8654      0 GAPDH   0     0   0   148.7583   148.7583   −149.7031      0.94421 PGK  1133.834     0  1133.834    0    0 2009588 −2016211 PGM   0      0   0    0    0   0     1959.931 EN   0      0   0    0    0    0      0 PK   2.067691     0   2.067691    0    0    0    −2.067691 LDH   0      0   0    0   0    39.25606      0 LEX   0      0   0    0    0    0      0 ATPase  −0.415409      0   −0.415409    0    0    0      0.415409 PGM EN PKLDH LEX ATPase HK     0   0     0.207471    0   0  −0.207471 PGI     0  0     0    0   0   0 PFK     0   0     0.207582    0   0  −0.207582ALD     0   0     0    0   0   0 TPI     0   0     0    0   0   0 GAPDH    0   0     0    0.00061   0   0 PGK    5490.018   0 −1133.834    0  0 1133.834 PGM −15287.46 13327.53     0    0   0   0 EN     81.36511 −81.44648     0.081365    0   0   0 PK     0   48.23862  −50.54352   0.237211   0   2.067691 LDH     0   0    232.001 −271.4036   0.146556  0 LEX     0   0     0    5.055766 −5.055766   0 ATPase     0   0    0.415409    0   0  −0.415409

The ability to describe the Jacobian matrix in terms of concentrationvariables or flux variables, and consequently being able to describe thedynamics in terms of fluxes or concentrations, is a result of theduality property described herein.

Example 4 Dynamic Properties of the Glycolytic Pathway in Homo sapiens

The construction of the key dynamic data matrices constructed for theglycolytic pathway in Homo sapiens in Example 1, 2, and 3, now enableanalysis of the dynamic properties of the pathway on different timescales.

TABLE 9 The metabolite modal matrix: Time Scales (hours) G6P F6P FDPDHAP GAP DPG13 PG3 4.95978E−07 −1.35273E−10   9.17E−09 −2.29E−09−1.04E−09   0.000127 −1.726916   0.004722  9.9213E−06  −0.411156665   1.003003  −3.7E−08   7.48E−13 −4.53E−09   3.06E−06 −8.54E−096.51335E−05 −4.34345E−08 −2.07E−08   4.25E−06   1.94E−06 −0.0018020.180719   0.17944   0.002226776  −0.001273406  −0.001254    0.105707   0.047837 −1.255173 −3.98E−05   8.46E−06   0.003681135  −0.00075682 −0.00075    0.026351    0.011856 −0.183504 −0.228332 −0.21689  0.019825246  −0.032068777  −0.032014    0.194885    0.072614 −0.149602−0.531172 −0.454169   0.027633412    0.381882846    0.381411  −2.204163   0.003881 −0.002006 −0.00655 −0.01758   0.104851825  −15.34226414 −15.33727   11.89888    6.163961 2.598787 −4.499937 −4.235859  0.141641663    70.63598759    70.61897 −22.40379 −11.47598 −6.562374  3.303189   2.598417   0.188045512 −267.7752441 −267.7267  −0.090144 −0.029436 −0.019904 −0.000565 −0.000398   0.197903124   215.6225745  215.5854   10.28596   5.212281   3.614971   0.386219   0.9089352668.439074    2.427898913    2.427899    3.236902   1.618453   1.618416  1.61834   0.809188 Time Scales (hours) PG2 PEP PYR LAC ATP 4.95978E−07−3.14E−05 −2.38E−08 −4.08E−10   1.77E−16   0.000974  9.9213E−06   1.3E−09 −1.39E−09   6.89E−12 −9.69E−18   2.91E−06 6.51335E−05−1.23007   6.87E−06   5.67E−09 −5.41E−14 −0.000102   0.002226776  8.67E−06   7.55E−06   4.24E−06  −1.4E−09 −4.77E−05   0.003681135−0.215398   0.059766 −1.426032 0.000784 −0.011439   0.019825246−0.442405   1.759867   0.001964 −6.36E−06 −0.07699   0.027633412−0.017255   0.043766   3.28E−05 −1.55E−07   0.01103   0.104851825−4.215199 −0.337184  −4.8E−05   1.62E−06 −0.26406   0.141641663  2.589042   0.828941    8.2E−05 −6.47E−06   0.704765   0.188045512−0.000397 −0.000188 −2.76E−08    3.5E−08 −0.000167   0.197903124 0.906966   0.537199   1.018456 0.999441 −0.522718 2668.439074  0.809188   0.809159   5.87E−09   1.65E−10   0.809152

TABLE 10 The flux modal matrix: Time Scales (hours) HK PGI PFK ALD TPIGAPDH PGK 4.95978E−07  −0.000974   9.31E−09  −0.000974   0.000127  0.000127 −1.727043   1.732612  9.9213E−06  −0.41116   1.41416 −1.003006   3.25E−08 −4.53E−09   3.06E−06 −1.53E−07 6.51335E−05   0.000102   2.28E−08    0.000106 −0.001805 −0.001804   0.182521−0.001381   0.002226776  −0.001226   1.95E−05    0.107008 −1.313044−1.30301   1.255134   5.99E−07   0.003681135    0.010682   7.02E−06   0.03854 −0.197999 −0.195361 −0.044828   2.97E−06   0.019825246   0.044921   5.52E−05    0.303888 −0.271873 −0.222216 −0.38157  1.32E−05   0.027633412    0.370853 −0.000472  −2.596604   2.206038−0.005888 −0.004544   1.16E−07   0.104851825  −15.0782   0.004993   27.50022 −3.136136 −3.565174 −7.098724    1.8E−05   0.141641663   69.93122 −0.017017  −93.72753   4.365443   4.913601   9.865563−6.96E−06   0.188045512 −267.7751   0.048591   267.6367   0.040804  0.009532   0.019339   1.06E−08   0.197903124   216.1453 −0.037179−204.7767 −1.45871 −1.59731 −3.228752 −2.12E−06 2668.439074    1.618747−2.55E−08  −0.000149 −3.36E−05 −3.68E−05 −7.56E−05 −3.37E−10 Time Scales(hours) PGM EN PK LDH LEX ATPase 4.95978E−07 −0.004754   3.14E−05  0.000974   4.08E−10 −1.77E−16 −0.000974  9.9213E−06   9.84E−09 −2.7E−09   2.91E−06 −6.89E−12   9.69E−18 −2.91E−06 6.51335E−05 −1.40951  1.230077 −0.000109 −5.67E−09   5.41E−14   0.000102   0.002226776  2.15E−07 −1.12E−06  −5.1E−05 −4.24E−06    1.4E−09   4.77E−05  0.003681135   0.001492   0.275164 −1.497237   1.426816 −0.000784  0.011439   0.019825246   0.011765   2.202271 −1.834892 −0.001971  6.36E−06   0.07699   0.027633412   0.000326   0.061021 −0.032703 −3.3E−05   1.55E−07 −0.01103   0.104851825   0.020659   3.878015  0.073076   4.96E−05 −1.62E−06   0.26406   0.141641663 −0.009375−1.760101 −0.124094 −8.85E−05   6.47E−06 −0.704765   0.188045512  1.11E−06   0.000209   2.13E−05   6.26E−08  −3.5E−08   0.000167  0.197903124 −0.001969 −0.369766 −0.041461 −0.019015 −0.999441  0.522718 2668.439074 −1.56E−07 −2.92E−05 −6.33E−06 −5.71E−09 −1.65E−10−0.809152

Pooling of variables (metabolites or fluxes) refers to the formation ofaggregate groups of variables, in which the group of variables movetogether in a concerted manner. Pools that form on very fast time scalesreflect the formation of chemical equilibrium pools, whereas poolingthat occurs on the slower time scales reflect physiologically relevantinteractions. Moving from the very fast time scales to the slower ones(FIG. 4), one observes the well-known examples of pool formation betweenhexose phosphates (HP) and phosphoglycerates (PG). With successiveremoval of modes on the slower time scales, more and more of themetabolites begin to form aggregate variables and move together in aconcerted fashion at fixed ratios. For glycolysis alone the successiveaggregation of chemical moieties (i.e. HP, PG) culminates in, on theslowest time scale, the formation of a physiologically meaningful poolthat represents the sum of high-energy phosphate bonds found in theglycolytic intermediates (i.e. their ATP equivalents) (FIGS. 3.4 and 4).The last row of M⁻¹ for J_(v) shows that this pool is moved byhexokinase as the input and ATPase as the output (FIG. 3B.4).

Example 5 Red Blood Cell Metabolism in Homo sapiens

The construction of the key dynamic data matrices constructed for a fullmodel of metabolites in Homo sapiens.

Currently, the human red cell is the only cell-scale kinetic modelavailable, whose formulation followed a 30 year history of iterativemodel building (summarized in (Joshi and Palsson, 1989; Joshi andPalsson, 1990), which are both hereby incorporated by reference in theirentireties). Analysis of the dynamic structure of this model resulted inthe simplification of network dynamics and the description of thecellular functions in terms of physiologically meaningful pools ofmetabolites. Analysis of the hierarchical dynamics of the human red cellmodel resulted in a richer and more complex physiological pool formation(Kauffman et al., 2002, which is hereby incorporated by reference) thatdetailed above for glycolysis alone. An overview of the hierarchicalreduction of the network into a functional diagram is schematized inFIG. 4. For example, the Adenosine Phosphate Potential is definedanalogously to Atkinson's Energy Charge (Atkinson, 1968, which is herebyincorporated by reference). As originally elaborated by Reich and Selkov(Reich and Selkov, 1981, which is hereby incorporated by reference),this “potential” is the ratio of the number of energy rich phosphatebonds and the ability or capacity to carry such bonds. The differentpools of metabolites in the red cell contribute to phosphate potentialsand/or oxidation/reduction potentials. The result of the pool formationis an ‘operating diagram’ (bottom of FIG. 5) that describes the functionof the metabolic network in the red cell slower time scales (minutes tohours).

This procedure characterized the dynamics of a network about aparticular steady state. The first step involves specifying a particularsteady state flux distribution. Next, the concentrations (generallythere will need to be some approximations of these) of the metabolitesat this steady state need to be identified. Third, the equilibriumconstant approximations need to be determined and the steady statefluxes and concentrations will be incorporated into mass action ratelaws. Finally, the forward rate constants can be solved by solving nlinear equations for forward rate constants.

For situations in which all of the concentrations are not known or ifone wants to trace the effects of particular metabolites in multiplerate constants, rather than n linear equations, one can solve the mequations by solving the steady state mass conservation relationship(S.v˜0). This will result in multiple alternative solutions.

As described previously (5), the dynamics of this model are consistentwith the observations from previous studies in which the well-knownexamples of pool formation between hexose phosphates (HP) andphosphoglycerates (PG). With successive removal of time scales (modes),more and more of the metabolites begin to form aggregate variables andmove together in a concerted fashion at fixed ratios. For glycolysisalone, the successive aggregation of chemical moieties (i.e. HP, PG)culminates in, on the slowest time scale, the formation of aphysiologically meaningful pool that represents the sum of high-energyphosphate bonds found in the glycolytic intermediates (i.e. their ATPequivalents). The last row of the modal matrix for the flux Jacobian(Supplemental Material) shows that this pool is moved by hexokinase asthe input and ATPase as the output, reflecting the catabolism of glucoseto generate ATP. The dimensions of the matrices and their correspondingsubspaces are shown in Table 11. This is a single pathway model, hencethe single dimension if the right null space is expected.

TABLE 11 The dimensions of the stoichiometric and gradient matrices,along with their ranks and sizes of their null spaces for the massaction glycolytic model. Rows Columns Rank Left Null Right NullStoichiometric Matrix 12 13 12 0 1 Gradient Matrix 13 12 12 0 0

Scaleability of the Approach

Since the approach outlined in FIG. 8 was successful for a small pathwaymodel we moved to apply this to expand the scope of the pathways.Building from the glycolytic pathway, the approach was applied to a fullmodel of human red blood cell metabolism. Steady state concentrationsfluxes and concentrations (16) reported in the literature were used inconjunction with equilibrium constants reported in the literature (17,18). The procedure outlined in FIG. 8 was again applied.

A map of the network is shown in FIG. 10 and the dimensions of thenetwork along with the size of the subspaces are shown in Table 12.

TABLE 12 The dimensions of the stoichiometric and gradient matrices,along with their ranks and sizes of their null spaces for the massaction model of the red cell network. Rows Columns Rank Left Null RightNull Stoichiometric Matrix 36 42 33 3 9 Gradient Matrix 42 36 36 6 0

As expected from the size of the left null space, there are threeconserved pools of metabolites, the stoichiometric sum of all of thephosphate containing compounds, the NAD moiety, and the NADP moiety. Theamount of glutathione in the network is also conserved, however sincethe dimeric form is not included as a dynamic variable in the network,it does not appear in the left null space. The size of the right nullspace increased significantly, by 8 dimensions, reflecting the variousalternate pathways that can occur from the Rappoport-Leubering Shunt,the pentose-phosphate pathway, and the connections between thenon-oxidative branches of the pentose phosphate pathway and thenucleotide salvage pathways.

The gradient and Jacobian matrices have the characteristics leadingtowards the expected behavior for the metabolic network, such ashierarchical structure with fast reactions corresponding toisomerization reactions and slower ones reflecting the morephysiological behaviors of the network, such as phosphate grouptransfers. There is good time scale separation (about 20 milliseconds to5 hours), although the slowest and fastest modes of the traditional redcell model is slightly larger (16, 19). These slower time scales in thekinetic model resulted from the potassium leak channels (20, 21), whichwere not included in this model. The compositions of the very fastmotions are similar to those seen for the red cell (e.g. isomerizationbetween glycolytic intermediates, hexose sugars, pentose sugars, etc).With this implementation, although the diphosphoglycerates move on theslower modes in the network, they are still faster than expected (on theorder of 5 minutes as opposed to 10 hours).

Construction of Mechanistic Regulated Model of the Human Erythrocyte

Fundamentally, all interactions between metabolites and enzymes arebilinear, including regulatory reactions (22). The mass action model forthe red blood cell was expanded in order in order to include anincreased level of detail for the biochemical reactions, accounting formechanistic regulation in the key regulatory enzymes of the red cell.Many of the association and dissociation constants for the stepsinvolved in enzymatic catalysis for the enzymes in the red cell havebeen characterized. The five key regulatory enzymes in the red bloodcell were described in terms of their mechanisms, these includehexokinase (HK), phosphofructokinase (PFK), diphosphoglyceromutate anddiphosphoglycerate phosphatase (DPGM and DPGase carried out by the sameenzyme), glucose-6-phosphate dehydrogenase (G6PDH), and adenosine kinase(AK). The mass action model discussed in the previous section was usedas a basis for this model. The incorporation of the bilinearinteractions reflecting small metabolite regulation of the enzymesresulted in a dramatic increase in the size of the network. Given thecomplexity and increase in size, the addition of each regulated enzymewas carried out in a systematic, ‘modular’ manner. Thus the parametersfor each enzyme were calculated individually. After each individualenzymatic step was completed and test, they were integrated into thewhole network model. Since a particular flux state was underconsideration, there were no problems with alternative fluxdistributions in the network.

The mechanism for unordered mechanism for hexokinase along with theassociation and dissociation constants were used as described in (23).FIG. 9 illustrates the individual steps. The total amount of enzyme isconserved (appears in the left null space) and was fixed at 24 nM asdescribed by (23). The net flux for the reaction was fixed according tothe steady state flux for the mass action model. The regulatoryinteractions by 2,3DPG and G6P do not carry net fluxes, so theirequilibrium concentrations were calculated. The remaining sets ofintermediate concentrations were then calculated.

PFK is a tetramer and as is apparent from the reaction diagram in FIG.9, it has four binding allosteric regulators. Although such interactionshave been described using the well known Hill equation, they can moreappropriately be directly described as a set of sequential bilinearinteractions. Association and dissociation constants for the regulatedsteps for PFK in the red blood cell were not available, hence anapproximation of 100,000 was made for most of the steps, for theexception of the catalytic transition, which was assumed to be muchslower. The total concentration of enzyme was 0.033 11M according to(24).

The diphosphoglyceromutase and diphosphoglycerol phosphatase reactionsare both carried out by the same enzyme and the enzymatic steps weredescribed according to (23). As with HK, the majority of the associationand dissociations steps were measured (23). A similar procedure wascarried out in order to calculate the concentrations of theenzyme-metabolite intermediates.

The reactions for G6PDH were implemented as described by (23) and solvedin an analogous manner as HK. The reaction scheme for AK is outlined inFIG. 9, as described by (25). The magnesium equilibrium reactions with2,3DPG and the adenosine phosphates were implemented as described in(17).

A summary of the dimensions and sizes of the subspaces for the regulatedmodel appear in Table 13.

TABLE 13 The dimensions of the stoichiometric and gradient matrices,along with their ranks and sizes of their null spaces for the regulatedmass action model of red cell metabolism. Rows Columns Rank Left NullRight Null Stoichiometric Matrix 92 94 82 10 12 Gradient Matrix 94 92 875 3

Note that this network is more than twice as large as the previousversion. Enzyme catalyzed mechanisms were introduced in this model for 5enzymes and the size of the left null space increased from 3 to 10. Thesize of the left null space grew by 7 dimensions due to the additionenzymes as reactants; 5 enzymes (note that PFK has an inactive andactive form) plus the conservation of total magnesium. The right nullspace increased by 3 as a result of the non-sequential mechanism of someof the regulated enzymes, such as HK.

Modal matrices for larger networks are much more difficult to interpretand were dependent on involved algorithmic approaches to elucidate thestructure (2). As alluded to earlier, this is because as the size ofnetworks increase and the non-integer entries in G which tends to be anillconditioned matrix, interactions become more tightly tied together.However, observing the structure of the modal matrix for the regulatednetwork (see Supplemental Material) a few points can be observed,

The time scale separation has improved from the non-regulated model,with motions ranging from less than microseconds to ten hours.

The slow modes are dominated by the metabolite bound enzymes.

In particular, the reactions that have a larger number of intermediatecomplexes playa more dominant role, especially at the slower timescales.

The movement of 2,3DPG (enzyme bound forms) have moved to the slowesttime scales.

In short, the structure of the modal matrix more closely approximatesthe features in the full red cell kinetic model. This has encouragingand interesting implications, since the starting point was a basic modelof glycolysis, expanded it to a model of red cell metabolism, and thenincorporated regulatory interactions in a mechanistic manner. Theresults from the regulated model of red cell metabolism suggest that onemanner in which regulatory enzymes exert their influence over the restof the network is their ability to bind more than one substrate and toexist in multiple states (active and inactive).

Simulations

Simulations were carried out with the regulated cell scale model inresponse to instantaneous ATP depletion (FIGS. 10 and 9) or NADPHdepletion. The large number of fluxes and concentrations in the networkpreclude the ability to plot all of them in a single figure, so keyfluxes and concentration variables in different pathways are plotted.Just as ratios of metabolites (e.g. Atkinson's Energy Charge) providemore insight than individual metabolite concentrations, we consider theratios of the enzymes in different states. In the simulations consideredin these studies, the ratio of free to unbound form of the regulatedenzymes. Under both simulation conditions, oscillatory motions betweenthe Rapport Leubering shunt, HK and PFK in glycolysis, and nucleotidemetabolism are observed. Based upon the large number of reactionsinvolved for the regulated steps as shown in FIG. 9 suggests that thereare a large number of possible interactions that can occur for eachindividual enzyme. This point is made explicitly clear in FIG. 11, whichdepicts the phase plane diagrams for all 21 of the reactions involvingphosphofructokinase. Some of the fluxes are highly correlated, whileothers exhibit independent motions along different time scales. Thus,regulated reactions are actually subnetworks within metabolic networks.

The NADHP/NADP ratio, G6PD enzyme ratio, and flux through the G6PDreaction are all highly correlated under both conditions (energy andredox stresses), suggesting that either of these measurements can beused to infer the states of the others. Correlations among many of theother variables appear to vary under different conditions.

Conclusions

Stoichiometric genome-scale models have made significant advancements inrecent years (26). The incorporation of dynamic information and theability to describe the kinetics of these networks is an important,physiologically relevant aspect that needs to be accounted for. Thekcone formalism took the initial steps towards developing kinetic modelsin a data driven fashion through sampling the kinetic parameter space(4). Herein we described an approach which uses stoichiometric models asa scaffold for constructing kinetic models based upon bilinear kinetics,through the incorporation of metabolomic data. This stepwise procedureis scaleable and feasible for networks of increasing size.

In order to further expand the mechanistic detail of stoichiometricmodels, we demonstrated how complex regulatory schemes, which aredescribed with more fidelity as a scheme of bilinear interactions, canfit within this framework of kinetic model construction. Subsequentanalysis of the network suggested that the regulatory influence ofenzymes is determined by the number of metabolites that can bind to theenzyme. We have previously shown that stoichiometry is a key aspect ofkinetic networks (5); here we further highlight that it is also a keyfactor in regulatory control as well.

Accounting for regulation allows one to observe the sub-networks ofregulation within regulated metabolic networks. These subnetworks canexert influences across the whole network, and this appears at least inpart in their ability to bind multiple compounds. As with smallmetabolites, the ratios of enzymes in different forms (e.g. bound versusunbound) reflect different functional states. For some enzymes these arehighly correlated with the flux and/or related small metabolite ratios,as with glucose-6-phosphate dehydrogenase. However for other enzymes,particularly those with a larger number of different binding states,such ratios will not always be correlated with fluxes.

Increasing availability of metabolomic data in conjunction with theavailability of Keq approximations and related methods for extrapolationand approximation of these constants (27-29) will enable the applicationof this approach to more metabolic networks. During the preparation ofthis work, Tran et al (10) describe an approach with similarities to thescheme described herein in which the solution space is sampled for aglycolysis/pentose phosphate pathway model in E. coli. The addition ofsampling would enable to the exploration of alternative solutions andhighlights an area worth further exploration in the future. Takentogether, the approach described within this work is scale-able andsystematic, thus representing a feasible, realistic path towardsgenome-scale kinetic descriptions of metabolic networks.

Methods

In general for any reaction, k,+

2A _(k1−←→) ^(k1+) B

the resulting net rate law is given by,

v]=k ₁ ⁺ A ₂ −k ₁ ⁻ B  (1),

with the equilibrium constant,

Keq=k ₁ ⁺ /k ₁ ⁻  (2)

For the regulated reactions in the network, the enzymes are explicitspecies and combined with each metabolite, as outlined in the figures.

Steps for Model Construction—Modular Model Construction

Substituting flux rate expressions (1) and equilibrium constants (2),into the steady state mass balances,

S·v=0  (3)

for which S is the m×n stoichiometric matrix and v is the flux vector oflength n, allows one to solve for the n forward rate constants k. Withthe rate constants and equilibrium constants in hand in addition tosteady state concentrations, the gradient matrix can be defined. Thegradient matrix G is given by dv/dx. Decomposition into the dualJacobian matrices (5) is given by

J _(x) −S·G=MΛM′;  (4)

m _(x) =M _(x) ⁻¹ x′  (5)

for the concentration Jacobian, J_(x), and

J _(v) =G·S=MΛM;′  (6)

m _(v) =M _(v) ⁻¹ v′  (7)

for the flux Jacobian, Jv. These calculations were carried out for eachof the sub-models that were developed in the studies described in themanuscript. The regulated cell model also included magnesium complexingwith ATP, ADP, AMP, and 2,3DPG. The matrices S and G are represented forvarious studies in FIGS. 13-16.

Models were constructed and tested in Mathematica (Wolfram Research).Simulations were also all carried out in Mathematica. However, modelconstruction and simulation can be carried out in any number of othercommercial software packages, such as Matlab, free software packagessuch as Scilab and R, or written in programming and scripting languagessuch as Java, C++, python, perl, etc. Since the protein concentrationswere much smaller than the concentrations of the metabolites,simulations were carried out using normalized enzyme concentrations, inorder to reduce the stiffness of the system of differential equations.

When carrying out simulations, the enzyme and bound enzyme variableswere normalized, by redefining the enzyme rate constant. This allowedthe stiffness of the ordinary differential equations to be significantlyreduced. The system was then integrated out to >1000 hours to define thesteady state.

This example shows how physiologically meaningful dynamic structuresform as a result of the particular numerical values in G and how theyoverlay on the network structure given by S. Not all sets of numericalvalues give this dynamic structure. It has been shown that geneticvariation, as represented by sequence polymorphrism in particularenzymes (pyruvate kinase and glucose-6-phosphate dehydrogenase), candisrupt this dynamic structure and lead to pathological states (Jamshidiet al., 2002).

Example 6 Simulating Genetic Deficiencies in Homo sapiens ThroughAltered Metabolite Measurements

Considering the dynamic model constructed in the above examples formetabolism in the human erythrocyte, one can simulate geneticdeficiencies through measurements of metabolites. For example a modelfor a cell with a deficiency in Phosphofructokinase can be constructedand simulated based on metabolite measurements. Following the algorithmfor model construction, a new set of forward rate constants will becalculated, resulting in a different set of reaction rate expressionsand gradient matrix. Subsequent simulations will reflect the resultingaltered biochemical phenotype due to the genetic mutation(s) resultingin Phosphofructokinase Deficiency.

Example 7 Simulating Genetic Deficiencies in Homo sapiens ThroughAltered Rate Constants

Considering the dynamic model constructed in the above examples formetabolism in the human erythrocyte, one can simulate geneticdeficiencies through measurements of kinetic parameters, such as rateconstants. For example, a model for a cell with a deficiency inPhosphofructokinase can be constructed and simulated based onmeasurements of altered kinetic parameters for the individual regulatedsteps. The mechanistic representation of the regulation of PFK by othermetabolites enables the representation of genetic mutations in a direct,detailed manner. Following the algorithm for model construction, thealtered rate constants are adjusted and an altered set of concentrationsare recalculated through solving the steady state equations andsimulating the time course of the dynamics. This will result in a newset of reaction rate expressions and gradient matrix. Subsequentsimulations will reflect the resulting altered biochemical phenotype dueto the genetic mutation(s) resulting in Phosphofructokinase Deficiency.

Example 8 Simulating Genetic Deficiencies in Homo sapiens ThroughAltered Lumped Parameter Measurements

Considering the dynamic model constructed in the above examples formetabolism in the human erythrocyte, one can simulate geneticdeficiencies through measurements of kinetic parameters, such as rateconstants. For example, a model for a cell with a deficiency inPhosphofructokinase can be constructed and simulated based onmeasurements of kinetic parameters such as V_(max). Following thealgorithm for model construction, the steady state equations and dynamicsimulations are re-implemented. This will result in a new set ofreaction rate expressions and gradient matrix. Subsequent simulationswill reflect the resulting altered biochemical phenotype due to thegenetic mutation(s) resulting in Phosphofructokinase Deficiency.

Example 9 Simulating the Effects of Drugs with Known Targets in Homosapiens

Considering the dynamic model constructed in the above examples formetabolism in the human erythrocyte, one can simulate the effects ofknown drugs on metabolism. For example, the effect of a redox depletingdrug which is known to act on a particular enzyme within the network,such as glucose-6-phosphate dehydrogenase, can be simulated by reducingor abrogating the rate constant for glucose-6-phosphate dehydrogenase by10%, 20%, 40%, 60%, 80%, or 100%. Following the algorithm for modelconstruction, the altered rate constants are adjusted and an altered setof concentrations are recalculated through solving the steady stateequations and simulating the time course of the dynamics. This willresult in a new set of reaction rate expressions and gradient matrix.The resulting model will reflect a metabolic network treated with thedrug at the specified dose.

Example 10 Simulating the Effects of Drugs with Unknown Targets in Homosapiens

Considering the dynamic model constructed in the above examples formetabolism in the human erythrocyte, one can simulate the effects ofknown drugs on metabolism and identify potential candidates for theidentity of the drug target. For example, metabolite concentrationmeasurements can be made experimentally following the administration ofa drug with unknown identity or unknown target. Following the algorithmfor model construction, a new set of forward rate constants will becalculated, resulting in a different set of reaction rate expressionsand gradient matrix. The resulting model will reflect a metabolicnetwork treated with the drug at the specified dose. Analysis of theJacobian matrices and dynamic simulations can yield insights into themechanism of action of the drug and its potential candidate targets.

REFERENCES

The following references are incorporated herein in their entirety:

-   1. Scriver, C. 2000. The Metabolic and Molecular Bases of Inherited    Disease. McGraw-Hill Professional.-   2. Jamshidi, N., and B. O. Palsson. 2008. Top-down analysis of    temporal hierarchy in biochemical reaction networks. PLoS Comput    BioI4:eI000177.-   3. Bruggeman, F. J., J. de Haan, H. Hardin, J. Bouwman, S.    Rossell, K. van Eunen, B. M. Bakker, and H. V. Westerhoff. 2006.    Time-dependent hierarchical regulation analysis: deciphering    cellular adaptation. Systems biology 153:318-322.-   4. Famili, I., R. Mahadevan, and B. O. Palsson. 2005. k-Cone    analysis: determining all candidate values for kinetic parameters on    a network scale. Biophys J 88:1616-1625.-   5. Jamshidi, N., and B. O. Palsson. 2008. Formulating genome-scale    kinetic models in the post-genome era. Molecular systems biology 4:    171.-   6. Maurya, M. R., and S. Subramaniam. 2007. A kinetic model for    calcium dynamics in RAW 264.7 cells: I. Mechanisms, parameters, and    subpopulational variability. Biophys J 93:709-728.-   7. Saucerman, J. J., L. L. Brunton, A. P. Michailova, and A. D.    McCulloch. 2003. Modeling beta-adrenergic control of cardiac myocyte    contractility in silico. The Journal of biological chemistry    278:47997-48003.-   8. Smallbone, K., E. Simeonidis, D. S. Broomhead, and D. B.    Kell. 2007. Something from nothing bridging the gap between    constraint-based and kinetic modelling FEBS J 274:5576-5585.-   9. Steuer, R., T. Gross, J. Selbig, and B. Blasius. 2006. Structural    kinetic modeling of metabolic networks. Proc Nat! Acad Sci USA 103:    11868-11873.-   10. Tran, L. M., M. L. Rizk, and J. C. Liao. 2008. Ensemble modeling    of metabolic networks. Biophys J 95:5606-5617.-   11. Wishart, D. S., D. Tzur, C. Knox, R. Eisner, A. C. Guo, N.    Young, D. Cheng, K. Jewell, D. Arndt, S. Sawhney, C. Fung, L.    Nikolai, M. Lewis, M. A. Coutouly, I. Forsythe, P. Tang, S.    Shrivastava, K. Jeroncic, P. Stothard, G. Amegbey, D. Block, D. D.    Hau, J. Wagner, J. Miniaci, M. Clements, M. Gebremedhin, N. Guo, Y.    Zhang, G. E. Duggan, G. D. Macinnis, A. M. Weljie, R.    Dowlatabadi, F. Bamforth, D. Clive, R. Greiner, L. Li, T.    Marrie, B. D. Sykes, H. J. Vogel, and L. Querengesser. 2007. HMDB:    the Human Metabolome Database. Nucleic acids research 35:D521-526.-   12. Bennett, B. D., J. Yuan, E. H. Kimball, and J. D.    Rabinowitz. 2008. Absolute quantitation of intracellular metabolite    concentrations by an isotope ratio-based approach. Nat Protoc    3:1299-1311-   13. Claus, S. P., T. M. Tsang, Y. Wang, O. Cloarec, E. Skordi, F. P.    Martin, S. Rezzi, A. Ross, S. Kochhar, E. Holmes, and J. K.    Nicholson. 2008. Systemic multicompartmental effects of the gut    microbiome on mouse metabolic phenotypes. Molecular systems biology    4:219.-   14. Ishii, N., K. Nakahigashi, T. Baba, M. Robert, T. Soga, A.    Kanai, T. Hirasawa, M. Naba, K. Hirai, A. Hogue, P. Y. Ho, Y.    Kakazu, K. Sugawara, S. Igarashi, S. Harada, T. Masuda, N.    Sugiyama, T. Togashi, M. Hasegawa, Y. Takai, K. Yugi, K. Arakawa, N.    Iwata, Y. Toya, Y. Nakayama, T. Nishioka, K. Shimizu, H. Mori,    and M. Tomita. 2007. Multiple high-throughput analyses monitor the    response of E. coli to perturbations. Science (New York, N.Y.    316:593-597.-   15. Smith, C. A., G. O'Maille, E. J. Want, C. Qin, S. A.    Trauger, T. R. Brandon, D. E. Custodio, R. Abagyan, and G.    Siuzdak. 2005. METLIN: a metabolite mass spectral database. Ther    Drug Moni t 27: 747-7 51.-   16. Joshi, A., and B. O. Palsson. 1990. Metabolic dynamics in the    human red cell. Part IV-Data prediction and some model computations.    J Theor Bioi 142:69-85.-   17. Joshi, A., and B. O. Palsson. 1990. Metabolic dynamics in the    human red cell. Part III—Metabolic reaction rates. J Theor Bioi    142:41-68.-   18. Reich, J., and E. Selkov. 1981. Energy metabolism of the cell: a    theoretical treatise. Academic Press.-   19. Jamshidi, N., J. S. Edwards, T. Fahland, G. M. Church, and B. O.    Palsson. 2001. Dynamic simulation of the human red blood cell    metabolic network. Bioinformatics 17:286-287.-   20. Jamshidi, N., and B. O. Palsson. 2006. Systems biology of the    human red blood cell. Blood cells, molecules & diseases 36:239-247.-   21. Kauffman, K. J., J. D. Pajerowski, N. Jamshidi, B. O. Palsson,    and J. S. Edwards. 2002. Description and analysis of metabolic    connectivity and dynamics in the human red blood cell. Biophys J    83:646-662.-   22. Segel, 1. 1975. Enzyme Kinetics. John Wiley & Sons.-   23. Mulquiney, P. J., and P. W. Kuchel. 1999. Model of    2,3-bisphosphoglycerate metabolism in the human erythrocyte based on    detailed enzyme kinetic equations: equations and parameter    refinement. Biochem J 342 Pt 3:581-596.-   24. Albe, K. R., M. H. Butler, and B. E. Wright. 1990. Cellular    concentrations of enzymes and their substrates. J Theor Bioi    143:163-195.-   25. Hawkins, C. F., and A. S. Bagnara. 1987. Adenosine kinase from    human erytluocytes: kinetic studies and characterization of    adenosine binding sites. Biochemistry 26: 1982-1987.-   26. Feist, A. M., M. J. Herrgard, 1. Thiele, J. L. Reed, and B. O.    Palsson. 2009. Reconstruction of biochemical networks in    microorganisms. Nat Rev Microbiol 7: 129-143.-   27. Alberty, R. A. 2003. Thermodynamics of biochemical reactions.    Wiley-Interscience, Hoboken, N.J.-   28. Jankowski, M. D., C. S. Henry, L. J. Broadbelt, and V.    Hatzimanikatis. 2008. Group contribution method for thermodynamic    analysis of complex metabolic networks. Biophys J 95:1487-1499.-   29. Mavrovouniotis, M. L. 1991. Estimation of standard Gibbs energy    changes of biotransformations. The Journal of biological chemistry    266:14440-14445.-   30. Mahadevan, R., and C. H. Schilling 2003. The effects of    alternate optimal solutions in constraint-based genome-scale    metabolic models. Metabolic engineering 5:264-276.-   31. Mulquiney, P. J., W. A. Bubb, and P. W. Kuchel. 1999. Model of    2,3-bisphosphoglycerate metabolism in the human erythrocyte based on    detailed enzyme kinetic equations: in vivo kinetic characterization    of 2,3-bisphosphoglycerate synthase/phosphatase using 13C and 31P    NMR. Biochem J 342 Pt 3:567-580.

What is claimed is:
 1. A computer-implemented method for developing adynamic network model of a biological network, the method comprising:providing, at a computer system, a network data structure relating aplurality of reactants to a plurality of reactions, the plurality ofreactions comprising one or more flux distributions, wherein each of thereactions comprises one or more reactants identified as a substrate ofthe reaction, one or more reactants identified as a product of thereaction and a stoichiometric coefficient relating the substrate and theproduct; providing, at a computer system, a constraint set for theplurality of reactions in said biological network associated with theplurality of reactions; providing, at a computer system, a constraintset for a plurality of constraints, wherein said constraint setcomprises at least one physical constraint; providing, at a computersystem, a constraint set for a plurality of kinetic constants associatedwith the plurality of reactions; providing, at a computer system, aconstraint set for a plurality of equilibrium constants associated withthe plurality of reactions; integrating, at a computer system, the setsof constraints to form stoichiometric, gradient, and Jacobian matrices;providing, at a computer system, an objective function; and predicting,at a computer system, a physiological function related to saidbiological network by determining using the computer system at least oneflux distribution, from the one or more flux distributions, thatminimizes or maximizes the objective function when at least oneconstraint set from the sets of constraints is applied to the datastructure.
 2. The method of claim 1, wherein said at least one physicalconstraint is mass balance.
 3. The method of claim 1, wherein saidconstraint set for said plurality of reactions in said biologicalnetwork reflects enzymatic interconversions of small metabolites, enzymereaction schemes, allosteric enzyme interactions, signaling pathways, ormacromolecular biochemical reactions.
 4. A computer-implemented methodfor developing a dynamic network model, the method comprising: accessingarchetype data for a plurality of reactions of a network; developing astoichiometric data structure for the reactions based on the archetypedata; accessing physiological data reflecting a specific context;performing a steady state simulation using the stoichiometric datastructure and the physiological data; incorporating metabolic data andequilibrium constant data based on the steady state simulation; derivinga gradient matrix and dynamic rate equations based on the stoichiometricdata structure and the incorporated data; and deriving a concentrationJacobian matrix and a flux Jacobian matrix based on the stoichiometricdata structure and the gradient matrix so as to characterize dynamicproperties of the network.
 5. The method of claim 4, wherein an i_(th)column of the stoichiometric data structure contains stoichiometriccoefficients of an i_(th) reaction in the network and an i_(th) row inthe gradient matrix contains partial derivatives of the i_(th) reactionwith respect to all reactants in the network.
 6. The method of claim 4,wherein the stoichiometric data structure is derived from genomic,proteomic and metabolomic data.
 7. The method of claim 4, wherein thegradient matrix is derived from fluxomic and metabolomic data, kineticcharacterization of individual reactions and assessment of thermodynamicproperties.
 8. The method of claim 4, wherein the stoichiometric datastructure is derived from a stoichiometry of reactions, and the gradientmatrix is based on kinetic data and thermodynamic information.
 9. Themethod of claim 4, wherein the stoichiometric data structure is based onthe content of a genome and is a property of a species, and the gradientmatrix is based on kinetic information and is genetically derived. 10.The method of claim 4, wherein a reaction comprises a representation ofa chemical conversion that consumes a substrate or produces a product,and a reactant comprises a representation of a chemical that is asubstrate or a product of a reaction that occurs in a cell.
 11. Themethod of claim 4, wherein a reaction comprises a conversion that occursdue to the activity of one or more enzymes that are genetically encodedby the human genome, a conversion that occurs spontaneously in a humanor mammalian cell, or conversions based on changes in chemicalcomposition such as those due to nucleophilic or electrophilic addition,nucleophilic or electrophilic substitution, elimination, isomerization,deamination, phosphorylation, methylation, glycolysation, reduction,oxidation or changes in location such as those that occur due to atransport reaction that moves one or more reactants within the samecompartment or from one cellular compartment to another.
 12. The methodof claim 4, wherein the biological network reflects a H. sapiens celltype at any stage of differentiation.
 13. The method of claim 4, whereinthe method is applied to normal cells or pathological cells.
 14. Themethod of claim 4, wherein the concentration Jacobian matrix iscalculated by post-multiplication of the stoichiometric data structureby the gradient matrix, and the flux Jacobian matrix is calculated bypre-multiplication of the stoichiometric data structure by the gradientmatrix.
 15. A computer-implemented method for developing data-drivendynamic models of biological networks, the method comprising: providing,at a computer system, a data structure relating a plurality ofbiological network reactants to a plurality of biological networkreactions, the plurality of biological network reactions comprising oneor more flux distributions, wherein each of the biological networkreactions comprises one or more biological network reactants identifiedas a substrate of the reaction, one or more biological network reactantsidentified as a product of the biological network reaction and astoichiometric coefficient relating the substrate and the product;providing, at a computer system, a constraint set for the plurality ofbiological network reactions; providing, at a computer system, aconstraint set for a plurality of concentrations associated with theplurality of biological network reactions; providing, at a computersystem, a constraint set for a plurality of kinetic constants associatedwith the plurality of biological network reactions; providing, at acomputer system, a constraint set for a plurality of equilibriumconstants associated with the plurality of biological network reactions;integrating, at a computer system, the sets of constraints, therebyforming stoichiometric, gradient, and Jacobian matrices; providing, at acomputer system, an objective function; and predicting, at a computersystem, a physiological function related to at least one of thebiological networks by determining using the computer system at leastone flux distribution, from the one or more flux distributions, thatminimizes or maximizes the objective function when at least oneconstraint set from the sets of constraints is applied to the datastructure, wherein one or more of said constraint sets reflects aperturbation.
 16. The method of claim 15, wherein one or more of saidconstraint sets has been altered to reflect genetic or environmentalperturbations.
 17. The method of claim 15, wherein at least one of saidkinetic constants or at least one of said concentrations has beenaltered to reflect a genetic perturbation.
 18. The method of claim 15,additionally comprising modifying, at a computer system, at least one ofthe reaction constraint set, the concentration constraint set, thekinetic constant set and the equilibrium constant set based at leastpartly on genetic variations.
 19. The method of claim 15, additionallycomprising modifying, at a computer system, at least one of the reactionconstraint set, the concentration constraint set, the constraint set fora plurality of kinetic constants and the set for a plurality ofequilibrium constants based at least partly on altered environmentalconditions.
 20. The method of claim 15, additionally comprisingsimulating the effects of adding or removing genes on at least one ofthe plurality of biological networks.
 21. The method of claim 15,additionally comprising simulating the effects of one or moreenvironmental perturbations on at least one of the plurality ofbiological networks.
 22. The method of claim 18, additionally comprisingdetermining dynamics and consequences of genetic variations selectedfrom the group consisting of deficiencies in metabolic enzymes andmetabolic transporters.
 23. The method of claim 22, wherein the geneticvariations comprise functioning of a compound selected from the groupconsisting of phosphofructokinase, phosphoglycerate kinase,phosphoglycerate mutase, lactate dehydrogenase adenosine deaminase, ABCtransporters, the SLC class of transporters, and the cytochrome P450class of enzymes.
 24. The method of claim 15, wherein providing aconstraint set for a plurality of concentrations comprises measuringmetabolite levels at steady state.
 25. The method of claim 15, whereinproviding a constraint set for a plurality of kinetic constantscomprises solving a steady state mass conservation relationshipcomprising a stoichiometric matrix and flux vector.
 26. The method ofclaim 19, wherein the modified sets of constraints correspond to anindividual.
 27. The method of claim 19, wherein the modified sets ofconstraints correspond to a group or plurality of individuals.
 28. Themethod of claim 18, wherein a therapeutic regime is applied to regulatea physiological function based on at least one of the modified setsbased at least partly on genetic variations.
 29. The method of claim 15,wherein a genetic variation results in alterations of kinetic parametersor concentrations of variables.
 30. The method of claim 19, wherein atherapeutic regime is applied to regulate a physiological function basedon at least one of the modified sets based at least partly on alteredenvironmental conditions.
 31. The method of claim 15, wherein thealtered environmental conditions result in changes in measuredvariables.
 32. The method of claim 15, additionally comprising:factorizing the gradient matrix into a kappa matrix and a gamma matrix;and analyzing the gradient, Jacobian, kappa and gamma matrices tocharacterize time scales of the network model.
 33. The method of claim32, wherein the kappa matrix is a diagonal matrix of kinetic constants.34. The method of claim 15, additionally comprising: modifying, at acomputer system, at least one of metabolite measurement or a kineticparameter measurement based on a genetic influence, the kineticparameter comprising a rate constant; generating a new set of reactionrate expressions and a new gradient matrix; and analyzing a pathologicalstate of the network model based on a simulation of the new set ofreaction rate expressions and the new gradient matrix.
 35. Acomputer-implemented method for the analysis of pathological states dueto genetic influences, the method comprising: providing, at a computersystem, a data structure relating a plurality of reactants to aplurality of reactions, the plurality of reactions comprising one ormore flux distributions, wherein each of the reactions comprises one ormore reactants identified as a substrate of the reaction, one or morereactants identified as a product of the reaction and a stoichiometriccoefficient relating the substrate and the product; providing, at acomputer system, a constraint set for the plurality of reactions;providing, at a computer system, a constraint set for a plurality ofconcentrations associated with the plurality of reactions; providing, ata computer system, a constraint set for a plurality of kinetic constantsassociated with the plurality of reactions; providing, at a computersystem, a constraint set for a plurality of equilibrium constantsassociated with the plurality of reactions; integrating, at a computersystem, the sets of constraints, thereby forming stoichiometric,gradient, and Jacobian matrices; providing, at a computer system, anobjective function; modifying, at a computer system, at least one of thereaction constraint set, the concentration constraint set, the kineticconstant set and the equilibrium constant set based at least partly ongenetic variations; and predicting, using a computer system, aphysiological function related to a pathological state by determining atleast one flux distribution, from the one or more flux distributions,that minimizes or maximizes the objective function when at least oneconstraint set from the sets of constraints is applied to the datastructure.
 36. The method of claim 35, additionally comprisingdetermining dynamics and consequences of genetic variations selectedfrom the group consisting of deficiencies in metabolic enzymes andmetabolic transporters.
 37. The method of claim 36, wherein the geneticvariations comprise functioning of a compound selected from the groupconsisting of phosphofructokinase, phosphoglycerate kinase,phosphoglycerate mutase, lactate dehydrogenase adenosine deaminase, ABCtransporters, the SLC class of transporters, and the cytochrome P450class of enzymes.
 38. The method of claim 35, wherein providing aconstraint set for a plurality of concentrations comprises measuringmetabolite levels at steady state.
 39. The method of claim 35, whereinproviding a constraint set for a plurality of kinetic constantscomprises solving a steady state mass conservation relationshipcomprising a stoichiometric matrix and flux vector.
 40. The method ofclaim 35, wherein the modified sets of constraints correspond to anindividual.
 41. The method of claim 35, wherein the modified sets ofconstraints correspond to a group or plurality of individuals.
 42. Themethod of claim 35, wherein a therapeutic regime is applied to regulatea physiological function based on at least one of the modified setsbased at least partly on genetic variations.
 43. The method of claim 35,wherein the genetic variations result in alterations of kineticparameters or concentrations of variables.
 44. A computer-implementedmethod for the analysis of pathological states due to environmentalinfluences, the method comprising: providing, at a computer system, adata structure relating a plurality of reactants to a plurality ofreactions, the plurality of reactions comprising one or more fluxdistributions, wherein each of the reactions comprises one or morereactants identified as a substrate of the reaction, one or morereactants identified as a product of the reaction and a stoichiometriccoefficient relating the substrate and the product; providing, at acomputer system, a constraint set for the plurality of reactions;providing, at a computer system, a constraint set for a plurality ofconcentrations associated with the plurality of reactions; providing, ata computer system, a constraint set for a plurality of kinetic constantsassociated with the plurality of reactions; providing, at a computersystem, a constraint set for a plurality of equilibrium constantsassociated with the plurality of reactions; integrating, at a computersystem, the sets of constraints, thereby forming stoichiometric,gradient, and Jacobian matrices; providing, at a computer system, anobjective function; modifying, at a computer system, at least one of thereaction constraint set, the concentration constraint set, theconstraint set for a plurality of kinetic constants and the set for aplurality of equilibrium constants based at least partly on alteredenvironmental conditions; and predicting, using a computer system, aphysiological function related to a pathological state by determining atleast one flux distribution, from the one or more flux distributions,that minimizes or maximizes the objective function when at least oneconstraint set from the sets of constraints is applied to the datastructure.
 45. The method of claim 44, wherein providing a constraintset for a plurality of concentrations comprises measuring metabolitelevels at steady state.
 46. The method of claim 44, wherein the modifiedsets of constraints correspond to an individual.
 47. The method of claim44, wherein the modified sets of constraints correspond to a group orplurality of individuals.
 48. The method of claim 44, wherein atherapeutic regime is applied to regulate a physiological function basedon at least one of the modified sets based at least partly on alteredenvironmental conditions.
 49. The method of claim 44, wherein thealtered environmental conditions result in changes in measuredvariables.
 50. A computer-implemented method for the analysis of apathological state of a data-driven model of a biological network due toa genetic influence, the method comprising: providing, at a computersystem, a data structure relating a plurality of reactants to aplurality of reactions, the plurality of reactions comprising one ormore flux distributions, wherein each of the reactions comprises one ormore reactants identified as a substrate of the reaction, one or morereactants identified as a product of the reaction and a stoichiometriccoefficient relating the substrate and the product; providing, at acomputer system, a constraint set for the plurality of reactions;providing, at a computer system, a constraint set for a plurality ofconcentrations associated with the plurality of reactions; providing, ata computer system, a constraint set for a plurality of kinetic constantsassociated with the plurality of reactions; providing, at a computersystem, a constraint set for a plurality of equilibrium constantsassociated with the plurality of reactions; integrating, at a computersystem, the sets of constraints, thereby forming stoichiometric,gradient, and Jacobian matrices; modifying, at a computer system, atleast one of metabolite measurement or a kinetic parameter measurementbased on a genetic influence, the kinetic parameter comprising a rateconstant; generating a new set of reaction rate expressions and a newgradient matrix; and analyzing of the pathological state of the networkmodel based on a simulation of the new set of reaction rate expressionsand the new gradient matrix.
 51. A computer-implemented method for theanalysis of an effect of a drug or chemical agent on a physiologicalfunction, the method comprising: providing, at a computer system, a datastructure relating a plurality of reactants to a plurality of reactions,the plurality of reactions comprising one or more flux distributions,wherein each of the reactions comprises one or more reactants identifiedas a substrate of the reaction, one or more reactants identified as aproduct of the reaction and a stoichiometric coefficient relating thesubstrate and the product; providing, at a computer system, a constraintset for the plurality of reactions; providing, at a computer system, aconstraint set for a plurality of concentrations associated with theplurality of reactions; providing, at a computer system, a constraintset for a plurality of kinetic constants associated with the pluralityof reactions; providing, at a computer system, a constraint set for aplurality of equilibrium constants associated with the plurality ofreactions; integrating, at a computer system, the sets of constraints,thereby forming the stoichiometric, gradient, and Jacobian matrices;providing, at a computer system, an objective function; adding areaction to the data structure, removing a reaction from the datastructure, or adjusting a constraint for a reaction in at least one ofthe sets of constraints to reflect a measured effect of a drug orchemical agent on the activity of the reaction; and predicting, at acomputer system, a physiological function related to the drug orchemical agent by determining at least one flux distribution, from theone or more flux distributions, that minimizes or maximizes theobjective function when a constraint set is applied to the datastructure.
 52. A computer-implemented method for developing a dynamicnetwork model, the method comprising: accessing archetype data for aplurality of reactions of a network; developing a stoichiometric datastructure for the reactions based on the archetype data; accessingphysiological data reflecting a specific context; performing a steadystate simulation using the stoichiometric data structure and thephysiological data; incorporating metabolic data and equilibriumconstant data based on the steady state simulation; deriving a gradientmatrix and dynamic rate equations based on the stoichiometric datastructure and the incorporated data; and deriving a concentrationJacobian matrix and a flux Jacobian matrix based on the stoichiometricdata structure and the gradient matrix so as to characterize dynamicproperties of the network, wherein at least some of the stoichiometric,metabolic or equilibrium constant data is obtained by sampling asolution space.
 53. A computer-implemented method for the analysis of aneffect of a drug or chemical agent on a physiological function, themethod comprising: accessing archetype data for a plurality of reactionsof a biological network; developing a stoichiometric data structure forthe reactions based on the archetype data; accessing physiological datareflecting a specific context; performing a steady state simulationusing the stoichiometric data structure and the physiological data;incorporating metabolic data and equilibrium constant data based on thesteady state simulation; adding a reaction to the stoichiometric datastructure, removing a reaction from the stoichiometric data structure,or adjusting a constraint for a reaction in at least a set ofconstraints for reactions in the stoichiometric data structure toreflect a measured effect of a drug or chemical agent on the activity ofthe reaction; deriving a gradient matrix and dynamic rate equationsbased on the stoichiometric data structure and the incorporated data;deriving a concentration Jacobian matrix and a flux Jacobian matrixbased on the stoichiometric data structure and the gradient matrix so asto characterize dynamic properties of the network; and predicting, at acomputer system, a physiological function related to the drug orchemical agent based on a simulation including at least the addedreaction, removed reaction, or adjusted constraint.
 54. Acomputer-implemented method for developing a dynamic network model, themethod comprising: accessing archetype data for a plurality of reactionsof a network; developing a stoichiometric data structure for thereactions based on the archetype data; performing a steady statesimulation using the stoichiometric data structure and physiologicaldata; obtaining metabolic data and equilibrium constant data based onthe steady state simulation; deriving a gradient matrix and dynamic rateequations based on the stoichiometric data structure and the obtaineddata; and deriving a concentration Jacobian matrix and a flux Jacobianmatrix based on the stoichiometric data structure and the gradientmatrix so as to characterize dynamic properties of the network, whereinunknown concentrations of metabolites are determined by solving a steadystate mass conservation relationship S·v=0, where S is a thestoichiometric data structure and v is a flux vector.